source/unres/io_config.f90

   1       module io_config
   2
   3       use names
   4       use io_units
   5       use io_base
   6       use geometry_data
   7       use geometry
   8       use control_data, only:maxterm_sccor
   9       implicit none
  10 !-----------------------------------------------------------------------------
  11 ! Max. number of residue types and parameters in expressions for
  12 ! virtual-bond angle bending potentials
  13 !      integer,parameter :: maxthetyp=3
  14 !      integer,parameter :: maxthetyp1=maxthetyp+1
  15 !     ,maxtheterm=20,
  16 !     & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
  17 !     & mmaxtheterm=maxtheterm)
  18 !-----------------------------------------------------------------------------
  19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
  20 ! and the number of terms in double torsionals
  21 !      integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
  22 !      parameter (maxtor=4,maxterm=10)
  23 !-----------------------------------------------------------------------------
  24 ! Max number of torsional terms in SCCOR
  25 !el      integer,parameter :: maxterm_sccor=6
  26 !-----------------------------------------------------------------------------
  27       character(len=1),dimension(:),allocatable :: secstruc     !(maxres)
  28 !-----------------------------------------------------------------------------
  29 !
  30 !
  31 !-----------------------------------------------------------------------------
  32       contains
  33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
  34 !-----------------------------------------------------------------------------
  35 ! bank.F    io_csa
  36 !-----------------------------------------------------------------------------
  37       subroutine write_rbank(jlee,adif,nft)
  38
  39       use csa_data
  40       use geometry_data, only: nres,rad2deg
  41 !      implicit real*8 (a-h,o-z)
  42 !      include 'DIMENSIONS'
  43 !      include 'COMMON.IOUNITS'
  44 !      include 'COMMON.CSA'
  45 !      include 'COMMON.BANK'
  46 !      include 'COMMON.CHAIN'
  47 !      include 'COMMON.GEO'
  48 !el local variables
  49       integer :: nft,i,k,j,l,jlee
  50       real(kind=8) :: adif
  51
  52       open(icsa_rbank,file=csa_rbank,status="unknown")
  53       write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
  54       do k=1,nbank
  55        write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
  56        do j=1,numch
  57         do l=2,nres-1
  58          write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
  59         enddo
  60        enddo
  61       enddo
  62       close(icsa_rbank)
  63
  64   850 format (10f8.3)
  65   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
  66               i8,i10,i2,f15.5)
  67   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
  68                   ' %NC ',0pf5.2)
  69
  70       return
  71       end subroutine write_rbank
  72 !-----------------------------------------------------------------------------
  73       subroutine read_rbank(jlee,adif)
  74
  75       use csa_data
  76       use geometry_data, only: nres,deg2rad
  77       use MPI_data
  78 !      implicit real*8 (a-h,o-z)
  79 !      include 'DIMENSIONS'
  80       include 'mpif.h'
  81 !      include 'COMMON.IOUNITS'
  82 !      include 'COMMON.CSA'
  83 !      include 'COMMON.BANK'
  84 !      include 'COMMON.CHAIN'
  85 !      include 'COMMON.GEO'
  86 !      include 'COMMON.SETUP'
  87       character(len=80) :: karta
  88 !el local variables
  89       integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
  90                  ierror,ierrcode,jlee,jleer
  91       real(kind=8) :: adif
  92
  93       open(icsa_rbank,file=csa_rbank,status="old")
  94       read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
  95       print *,jleer,nbankr,nstepr,nftr,icycler,adif
  96 !       print *, 'adif from read_rbank ',adif
  97 #ifdef MPI
  98       if(nbankr.ne.nbank) then
  99         write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
 100         ' NBANK',nbank
 101         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 102       endif
 103       if(jleer.ne.jlee) then
 104         write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
 105         ' JLEE',jlee
 106         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 107       endif
 108 #endif
 109
 110       kk=0
 111       do k=1,nbankr
 112         read (icsa_rbank,'(a80)') karta
 113         write(iout,*) "READ_RBANK: kk=",kk
 114         write(iout,*) karta
 115 !        if (index(karta,"*").gt.0) then
 116 !          write (iout,*) "***** Stars in bankr ***** k=",k,
 117 !     &      " skipped"
 118 !          do j=1,numch
 119 !            do l=2,nres-1
 120 !              read (30,850) (rdummy,i=1,4)
 121 !            enddo
 122 !          enddo
 123 !        else
 124           kk=kk+1
 125           call reada(karta,"total E",rene(kk),1.0d20)
 126           call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
 127           call reada(karta,"%NC",rpncn(kk),0.0d0)
 128           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
 129             "%NC",bpncn(kk),ibank(kk)
 130 !          read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
 131           do j=1,numch
 132             do l=2,nres-1
 133               read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
 134 !              write (iout,850) (rvar(i,l,j,kk),i=1,4)
 135               do i=1,4
 136                 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
 137               enddo
 138             enddo
 139           enddo
 140 !        endif
 141       enddo
 142 !d      write (*,*) "read_rbank ******************* kk",kk,
 143 !d     &  "nbankr",nbankr
 144       if (kk.lt.nbankr) nbankr=kk
 145 !d      do kk=1,nbankr
 146 !d          print *,"kk=",kk
 147 !d          do j=1,numch
 148 !d            do l=2,nres-1
 149 !d              write (*,850) (rvar(i,l,j,kk),i=1,4)
 150 !d            enddo
 151 !d          enddo
 152 !d      enddo
 153       close(icsa_rbank)
 154
 155   850 format (10f8.3)
 156   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 157   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
 158
 159       return
 160       end subroutine read_rbank
 161 !-----------------------------------------------------------------------------
 162       subroutine write_bank(jlee,nft)
 163
 164       use csa_data
 165       use control_data, only: vdisulf
 166       use geometry_data, only: nres,rad2deg
 167 !      implicit real*8 (a-h,o-z)
 168 !      include 'DIMENSIONS'
 169 !      include 'COMMON.IOUNITS'
 170 !      include 'COMMON.CSA'
 171 !      include 'COMMON.BANK'
 172 !      include 'COMMON.CHAIN'
 173 !      include 'COMMON.GEO'
 174 !      include 'COMMON.SBRIDGE'
 175 !     include 'COMMON.CONTROL'
 176       character(len=7) :: chtmp
 177       character(len=40) :: chfrm
 178 !el      external ilen
 179 !el local variables
 180       integer :: nft,k,l,i,j,jlee
 181
 182       open(icsa_bank,file=csa_bank,status="unknown")
 183       write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
 184       write (icsa_bank,902) nglob_csa, eglob_csa
 185       open (igeom,file=intname,status='UNKNOWN')
 186       do k=1,nbank
 187        write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 188        if (vdisulf) write (icsa_bank,'(101i4)') &
 189           bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
 190        do j=1,numch
 191         do l=2,nres-1
 192          write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 193         enddo
 194        enddo
 195        if (bvar_nss(k).le.9) then
 196          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 197            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 198        else
 199          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 200            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 201          write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
 202                                       bvar_ss(2,i,k),i=10,bvar_nss(k))
 203        endif
 204        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 205        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 206        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 207        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 208       enddo
 209       close(icsa_bank)
 210       close(igeom)
 211
 212       if (nstep/200.gt.ilastnstep) then
 213
 214        ilastnstep=(ilastnstep+1)*1.5
 215        write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
 216        write(chtmp,chfrm) nstep
 217        open(icsa_int,file=prefix(:ilen(prefix)) &
 218                //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
 219        do k=1,nbank
 220         if (bvar_nss(k).le.9) then
 221          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 222         bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 223         else
 224          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 225            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 226          write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
 227                                 bvar_ss(2,i,k),i=10,bvar_nss(k))
 228         endif
 229         write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 230         write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 231         write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 232         write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 233        enddo
 234        close(icsa_int)
 235       endif
 236
 237
 238   200 format (8f10.4)
 239   850 format (10f8.3)
 240   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
 241               i8,i10,i2,f15.5)
 242   902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
 243   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 244               ' %NC ',0pf5.2,i5)
 245
 246       return
 247       end subroutine write_bank
 248 !-----------------------------------------------------------------------------
 249       subroutine write_bank_reminimized(jlee,nft)
 250
 251       use csa_data
 252       use geometry_data, only: nres,rad2deg
 253       use energy_data
 254 !      implicit real*8 (a-h,o-z)
 255 !      include 'DIMENSIONS'
 256 !      include 'COMMON.IOUNITS'
 257 !      include 'COMMON.CSA'
 258 !      include 'COMMON.BANK'
 259 !      include 'COMMON.CHAIN'
 260 !      include 'COMMON.GEO'
 261 !      include 'COMMON.SBRIDGE'
 262 !el local variables
 263       integer :: nft,i,l,j,k,jlee
 264
 265       open(icsa_bank_reminimized,file=csa_bank_reminimized,&
 266         status="unknown")
 267       write (icsa_bank_reminimized,900) &
 268         jlee,nbank,nstep,nft,icycle,cutdif
 269       open (igeom,file=intname,status='UNKNOWN')
 270       do k=1,nbank
 271        write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
 272         bpncn(k),ibank(k)
 273        do j=1,numch
 274         do l=2,nres-1
 275          write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 276         enddo
 277        enddo
 278        if (nss.le.9) then
 279          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 280            nss,(ihpb(i),jhpb(i),i=1,nss)
 281        else
 282          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 283            nss,(ihpb(i),jhpb(i),i=1,9)
 284          write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
 285        endif
 286        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 287        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 288        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 289        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 290       enddo
 291       close(icsa_bank_reminimized)
 292       close(igeom)
 293
 294   200 format (8f10.4)
 295   850 format (10f8.3)
 296   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
 297               i8,i10,i2,f15.5)
 298   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 299                ' %NC ',0pf5.2,i5)
 300
 301       return
 302       end subroutine write_bank_reminimized
 303 !-----------------------------------------------------------------------------
 304       subroutine read_bank(jlee,nft,cutdifr)
 305
 306       use csa_data
 307       use control_data, only: vdisulf
 308       use geometry_data, only: nres,deg2rad
 309       use energy_data
 310 !      implicit real*8 (a-h,o-z)
 311 !      include 'DIMENSIONS'
 312 !      include 'COMMON.IOUNITS'
 313 !      include 'COMMON.CSA'
 314 !      include 'COMMON.BANK'
 315 !      include 'COMMON.CHAIN'
 316 !      include 'COMMON.GEO'
 317 !      include 'COMMON.CONTROL'
 318 !      include 'COMMON.SBRIDGE'
 319       character(len=80) :: karta
 320 !      integer ilen
 321 !el      external ilen
 322 !el local variables
 323       integer :: nft,kk,k,l,i,j,jlee
 324       real(kind=8) :: cutdifr
 325
 326       open(icsa_bank,file=csa_bank,status="old")
 327        read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
 328        read (icsa_bank,902) nglob_csa, eglob_csa
 329 !      if(jleer.ne.jlee) then
 330 !        write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
 331 !    &   ' JLEE',jlee
 332 !        call mpi_abort(mpi_comm_world,ierror,ierrcode)
 333 !      endif
 334
 335       kk=0
 336       do k=1,nbank
 337         read (icsa_bank,'(a80)') karta
 338         write(iout,*) "READ_BANK: kk=",kk
 339         write(iout,*) karta
 340 !        if (index(karta,"*").gt.0) then
 341 !          write (iout,*) "***** Stars in bank ***** k=",k,
 342 !     &      " skipped"
 343 !          do j=1,numch
 344 !            do l=2,nres-1
 345 !              read (33,850) (rdummy,i=1,4)
 346 !            enddo
 347 !          enddo
 348 !        else
 349           kk=kk+1
 350           call reada(karta,"total E",bene(kk),1.0d20)
 351           call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
 352           call reada(karta,"%NC",bpncn(kk),0.0d0)
 353           read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
 354           goto 112
 355   111     ibank(kk)=0
 356   112     continue
 357           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
 358             "%NC",bpncn(kk),ibank(kk)
 359 !          read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 360           if (vdisulf) then
 361             read (icsa_bank,'(101i4)') &
 362               bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 363             bvar_ns(kk)=ns-2*bvar_nss(kk)
 364             write(iout,*) 'read SSBOND',bvar_nss(kk),&
 365                           ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 366 !d          write(iout,*) 'read CYS #free ', bvar_ns(kk)
 367             l=0
 368             do i=1,ns
 369              j=1
 370              do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
 371                        iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
 372                        j.le.bvar_nss(kk))
 373                 j=j+1
 374              enddo
 375              if (j.gt.bvar_nss(kk)) then
 376                l=l+1
 377                bvar_s(l,kk)=iss(i)
 378              endif
 379             enddo
 380 !d            write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
 381           endif
 382           do j=1,numch
 383             do l=2,nres-1
 384               read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
 385 !              write (iout,850) (bvar(i,l,j,kk),i=1,4)
 386               do i=1,4
 387                 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
 388               enddo ! l
 389             enddo ! l
 390           enddo ! j
 391 !        endif
 392       enddo ! k
 393
 394       if (kk.lt.nbank) nbank=kk
 395 !d      write (*,*) "read_bank ******************* kk",kk,
 396 !d     &  "nbank",nbank
 397 !d      do kk=1,nbank
 398 !d          print *,"kk=",kk
 399 !d          do j=1,numch
 400 !d            do l=2,nres-1
 401 !d              write (*,850) (bvar(i,l,j,kk),i=1,4)
 402 !d            enddo
 403 !d          enddo
 404 !d      enddo
 405
 406 !       do k=1,nbank
 407 !        read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 408 !        do j=1,numch
 409 !         do l=2,nres-1
 410 !          read (33,850) (bvar(i,l,j,k),i=1,4)
 411 !          do i=1,4
 412 !           bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
 413 !          enddo
 414 !         enddo
 415 !        enddo
 416 !       enddo
 417       close(icsa_bank)
 418
 419   850 format (10f8.3)
 420   952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
 421   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 422   902 format (1x,11x,i4,12x,1pe14.5)
 423   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
 424
 425       return
 426       end subroutine read_bank
 427 !-----------------------------------------------------------------------------
 428       subroutine write_bank1(jlee)
 429
 430       use csa_data
 431       use geometry_data, only: nres,rad2deg
 432 !      implicit real*8 (a-h,o-z)
 433 !      include 'DIMENSIONS'
 434 !      include 'COMMON.IOUNITS'
 435 !      include 'COMMON.CSA'
 436 !      include 'COMMON.BANK'
 437 !      include 'COMMON.CHAIN'
 438 !      include 'COMMON.GEO'
 439 !el local variables
 440       integer :: k,i,l,j,jlee
 441
 442 #if defined(AIX) || defined(PGI)
 443       open(icsa_bank1,file=csa_bank1,position="append")
 444 #else
 445       open(icsa_bank1,file=csa_bank1,access="append")
 446 #endif
 447       write (icsa_bank1,900) jlee,nbank,nstep,cutdif
 448       do k=1,nbank
 449        write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 450        do j=1,numch
 451         do l=2,nres-1
 452          write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 453         enddo
 454        enddo
 455       enddo
 456       close(icsa_bank1)
 457   850 format (10f8.3)
 458   900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
 459   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 460               ' %NC ',0pf5.2,i5)
 461
 462       return
 463       end subroutine write_bank1
 464 !-----------------------------------------------------------------------------
 465 ! cartprint.f
 466 !-----------------------------------------------------------------------------
 467 !      subroutine cartprint
 468
 469 !      use geometry_data, only: c
 470 !      use energy_data, only: itype
 471 !      implicit real*8 (a-h,o-z)
 472 !      include 'DIMENSIONS'
 473 !      include 'COMMON.CHAIN'
 474 !      include 'COMMON.INTERACT'
 475 !      include 'COMMON.NAMES'
 476 !      include 'COMMON.IOUNITS'
 477 !      integer :: i
 478
 479 !     write (iout,100)
 480 !      do i=1,nres
 481 !        write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
 482 !          c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
 483 !      enddo
 484 !  100 format (//'              alpha-carbon coordinates       ',&
 485 !                '     centroid coordinates'/ &
 486 !                '       ', 6X,'X',11X,'Y',11X,'Z',&
 487 !                                10X,'X',11X,'Y',11X,'Z')
 488 !  110 format (a,'(',i3,')',6f12.5)
 489 !      return
 490 !      end subroutine cartprint
 491 !-----------------------------------------------------------------------------
 492 ! dihed_cons.F
 493 !-----------------------------------------------------------------------------
 494       subroutine secstrp2dihc
 495
 496       use geometry_data
 497       use energy_data
 498 !      implicit real*8 (a-h,o-z)
 499 !      include 'DIMENSIONS'
 500 !      include 'COMMON.GEO'
 501 !      include 'COMMON.BOUNDS'
 502 !      include 'COMMON.CHAIN'
 503 !      include 'COMMON.TORCNSTR'
 504 !      include 'COMMON.IOUNITS'
 505 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 506 !el      COMMON/SECONDARYS/secstruc
 507       character(len=80) :: line
 508       logical :: errflag
 509 !el      external ilen
 510
 511 !el  local variables
 512       integer :: i,ii,lenpre
 513
 514       allocate(secstruc(nres))
 515
 516 !dr      call getenv_loc('SECPREDFIL',secpred)
 517       lenpre=ilen(prefix)
 518       secpred=prefix(:lenpre)//'.spred'
 519
 520 #if defined(WINIFL) || defined(WINPGI)
 521       open(isecpred,file=secpred,status='old',readonly,shared)
 522 #elif (defined CRAY) || (defined AIX)
 523       open(isecpred,file=secpred,status='old',action='read')
 524 #elif (defined G77)
 525       open(isecpred,file=secpred,status='old')
 526 #else
 527       open(isecpred,file=secpred,status='old',action='read')
 528 #endif
 529 ! read secondary structure prediction from JPRED here!
 530 !      read(isecpred,'(A80)',err=100,end=100) line
 531 !      read(line,'(f10.3)',err=110) ftors
 532        read(isecpred,'(f10.3)',err=110) ftors
 533
 534       write (iout,*) 'FTORS factor =',ftors
 535 ! initialize secstruc to any
 536        do i=1,nres
 537         secstruc(i) ='-'
 538       enddo
 539       ndih_constr=0
 540       ndih_nconstr=0
 541
 542       call read_secstr_pred(isecpred,iout,errflag)
 543       if (errflag) then
 544          write(iout,*)'There is a problem with the list of secondary-',&
 545            'structure prediction'
 546          goto 100
 547       endif
 548 ! 8/13/98 Set limits to generating the dihedral angles
 549       do i=1,nres
 550         phibound(1,i)=-pi
 551         phibound(2,i)=pi
 552       enddo
 553
 554       ii=0
 555       do i=1,nres
 556         if ( secstruc(i) .eq. 'H') then
 557 ! Helix restraints for this residue
 558            ii=ii+1
 559            idih_constr(ii)=i
 560            phi0(ii) = 45.0D0*deg2rad
 561            drange(ii)= 5.0D0*deg2rad
 562            phibound(1,i) = phi0(ii)-drange(ii)
 563            phibound(2,i) = phi0(ii)+drange(ii)
 564         else if (secstruc(i) .eq. 'E') then
 565 ! strand restraints for this residue
 566            ii=ii+1
 567            idih_constr(ii)=i
 568            phi0(ii) = 180.0D0*deg2rad
 569            drange(ii)= 5.0D0*deg2rad
 570            phibound(1,i) = phi0(ii)-drange(ii)
 571            phibound(2,i) = phi0(ii)+drange(ii)
 572         else
 573 ! no restraints for this residue
 574            ndih_nconstr=ndih_nconstr+1
 575            idih_nconstr(ndih_nconstr)=i
 576         endif
 577       enddo
 578       ndih_constr=ii
 579 !      deallocate(secstruc)
 580       return
 581 100   continue
 582       write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
 583 !      deallocate(secstruc)
 584       return
 585 110   continue
 586       write(iout,'(A20)')'Error reading FTORS'
 587 !      deallocate(secstruc)
 588       return
 589       end subroutine secstrp2dihc
 590 !-----------------------------------------------------------------------------
 591       subroutine read_secstr_pred(jin,jout,errors)
 592
 593 !      implicit real*8 (a-h,o-z)
 594 !      INCLUDE 'DIMENSIONS'
 595 !      include 'COMMON.IOUNITS'
 596 !      include 'COMMON.CHAIN'
 597 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 598 !el      COMMON/SECONDARYS/secstruc
 599 !el      EXTERNAL ILEN
 600       character(len=80) :: line,line1   !,ucase
 601       logical :: errflag,errors,blankline
 602
 603 !el  local variables
 604       integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
 605             length_of_chain
 606       errors=.false.
 607       read (jin,'(a)') line
 608       write (jout,'(2a)') '> ',line(1:78)
 609       line1=ucase(line)
 610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
 611 ! correspond to the end-groups.  ADD to the secondary structure prediction "-" for the
 612 ! end-groups in the input file "*.spred"
 613
 614       iseq=1
 615       do while (index(line1,'$END').eq.0)
 616 ! Override commented lines.
 617          ipos=1
 618          blankline=.false.
 619          do while (.not.blankline)
 620             line1=' '
 621             call mykey(line,line1,ipos,blankline,errflag)
 622             if (errflag) write (jout,'(2a)') &
 623        'Error when reading sequence in line: ',line
 624             errors=errors .or. errflag
 625             if (.not. blankline .and. .not. errflag) then
 626                ipos1=2
 627                iend=ilen(line1)
 628 !el               if (iseq.le.maxres) then
 629                   if (line1(1:1).eq.'-' ) then
 630                      secstruc(iseq)=line1(1:1)
 631                   else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
 632                             ( ucase(line1(1:1)).eq.'H' ) ) then
 633                      secstruc(iseq)=ucase(line1(1:1))
 634                   else
 635                      errors=.true.
 636                      write (jout,1010) line1(1:1), iseq
 637                      goto 80
 638                   endif
 639 !el               else
 640 !el                  errors=.true.
 641 !el                  write (jout,1000) iseq,maxres
 642 !el                  goto 80
 643 !el               endif
 644                do while (ipos1.le.iend)
 645
 646                   iseq=iseq+1
 647                   il=1
 648                   ipos1=ipos1+1
 649 !el                  if (iseq.le.maxres) then
 650                      if (line1(ipos1-1:ipos1-1).eq.'-' ) then
 651                         secstruc(iseq)=line1(ipos1-1:ipos1-1)
 652                      else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
 653                            (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
 654                         secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
 655                      else
 656                         errors=.true.
 657                         write (jout,1010) line1(ipos1-1:ipos1-1), iseq
 658                         goto 80
 659                      endif
 660 !el                  else
 661 !el                     errors=.true.
 662 !el                     write (jout,1000) iseq,maxres
 663 !el                     goto 80
 664 !el                  endif
 665                enddo
 666                iseq=iseq+1
 667             endif
 668          enddo
 669          read (jin,'(a)') line
 670          write (jout,'(2a)') '> ',line(1:78)
 671          line1=ucase(line)
 672       enddo
 673
 674 !d    write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
 675
 676 !d check whether the found length of the chain is correct.
 677       length_of_chain=iseq-1
 678       if (length_of_chain .ne. nres) then
 679 !        errors=.true.
 680         write (jout,'(a,i4,a,i4,a)') &
 681        'Error: the number of labels specified in $SEC_STRUC_PRED (' &
 682        ,length_of_chain,') does not match with the number of residues (' &
 683        ,nres,').'
 684       endif
 685    80 continue
 686
 687  1000 format('Error - the number of residues (',i4,&
 688        ') has exceeded maximum (',i4,').')
 689  1010 format ('Error - unrecognized secondary structure label',a4,&
 690        ' in position',i4)
 691       return
 692       end subroutine read_secstr_pred
 693 !#endif
 694 !-----------------------------------------------------------------------------
 695 ! parmread.F
 696 !-----------------------------------------------------------------------------
 697       subroutine parmread
 698
 699       use geometry_data
 700       use energy_data
 701       use control_data, only:maxterm !,maxtor
 702       use MD_data
 703       use MPI_data
 704 !el      use map_data
 705       use control, only: getenv_loc
 706 !
 707 ! Read the parameters of the probability distributions of the virtual-bond
 708 ! valence angles and the side chains and energy parameters.
 709 !
 710 ! Important! Energy-term weights ARE NOT read here; they are read from the
 711 ! main input file instead, because NO defaults have yet been set for these
 712 ! parameters.
 713 !
 714 !      implicit real*8 (a-h,o-z)
 715 !      include 'DIMENSIONS'
 716 #ifdef MPI
 717       include "mpif.h"
 718       integer :: IERROR
 719 #endif
 720 !      include 'COMMON.IOUNITS'
 721 !      include 'COMMON.CHAIN'
 722 !      include 'COMMON.INTERACT'
 723 !      include 'COMMON.GEO'
 724 !      include 'COMMON.LOCAL'
 725 !      include 'COMMON.TORSION'
 726 !      include 'COMMON.SCCOR'
 727 !      include 'COMMON.SCROT'
 728 !      include 'COMMON.FFIELD'
 729 !      include 'COMMON.NAMES'
 730 !      include 'COMMON.SBRIDGE'
 731 !      include 'COMMON.MD'
 732 !      include 'COMMON.SETUP'
 733       character(len=1) :: t1,t2,t3
 734       character(len=1) :: onelett(4) = (/"G","A","P","D"/)
 735       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
 736       logical :: lprint,LaTeX
 737       real(kind=8),dimension(3,3,maxlob) :: blower      !(3,3,maxlob)
 738       real(kind=8),dimension(13) :: b
 739       character(len=3) :: lancuch       !,ucase
 740 !el  local variables
 741       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
 742       integer :: maxinter,junk,kk,ii,ncatprotparm
 743       real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
 744                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
 745                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
 746                 res1,epsijlip,epspeptube,epssctube,sigmapeptube,      &
 747                 sigmasctube
 748       integer :: ichir1,ichir2
 749 !      real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
 750 !el      allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
 751 !el      allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
 752 !
 753 ! For printing parameters after they are read set the following in the UNRES
 754 ! C-shell script:
 755 !
 756 ! setenv PRINT_PARM YES
 757 !
 758 ! To print parameters in LaTeX format rather than as ASCII tables:
 759 !
 760 ! setenv LATEX YES
 761 !
 762       call getenv_loc("PRINT_PARM",lancuch)
 763       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
 764       call getenv_loc("LATEX",lancuch)
 765       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
 766 !
 767       dwa16=2.0d0**(1.0d0/6.0d0)
 768       itypro=20
 769 ! Assign virtual-bond length
 770       vbl=3.8D0
 771       vblinv=1.0D0/vbl
 772       vblinv2=vblinv*vblinv
 773 !
 774 ! Read the virtual-bond parameters, masses, and moments of inertia
 775 ! and Stokes' radii of the peptide group and side chains
 776 !
 777       allocate(dsc(ntyp1)) !(ntyp1)
 778       allocate(dsc_inv(ntyp1)) !(ntyp1)
 779       allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
 780       allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
 781       allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
 782       allocate(nbondterm(ntyp)) !(ntyp)
 783       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 784       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 785       allocate(msc(ntyp+1,5)) !(ntyp+1)
 786       allocate(isc(ntyp+1,5)) !(ntyp+1)
 787       allocate(restok(ntyp+1,5)) !(ntyp+1)
 788       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 789       allocate(long_r_sidechain(ntyp))
 790       allocate(short_r_sidechain(ntyp))
 791       dsc(:)=0.0d0
 792       dsc_inv(:)=0.0d0
 793
 794 #ifdef CRYST_BOND
 795       read (ibond,*) vbldp0,akp,mp,ip,pstok
 796       do i=1,ntyp
 797         nbondterm(i)=1
 798         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
 799         dsc(i) = vbldsc0(1,i)
 800         if (i.eq.10) then
 801           dsc_inv(i)=0.0D0
 802         else
 803           dsc_inv(i)=1.0D0/dsc(i)
 804         endif
 805       enddo
 806 #else
 807       read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
 808       do i=1,ntyp_molec(1)
 809         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
 810          j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
 811         dsc(i) = vbldsc0(1,i)
 812         if (i.eq.10) then
 813           dsc_inv(i)=0.0D0
 814         else
 815           dsc_inv(i)=1.0D0/dsc(i)
 816         endif
 817       enddo
 818 #endif
 819       read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
 820       do i=1,ntyp_molec(2)
 821         nbondterm_nucl(i)=1
 822         read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
 823 !        dsc(i) = vbldsc0_nucl(1,i)
 824 !        if (i.eq.10) then
 825 !          dsc_inv(i)=0.0D0
 826 !        else
 827 !          dsc_inv(i)=1.0D0/dsc(i)
 828 !        endif
 829       enddo
 830 !      read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
 831 !      do i=1,ntyp_molec(2)
 832 !        read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
 833 !        aksc_nucl(j,i),abond0_nucl(j,i),&
 834 !         j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
 835 !        dsc(i) = vbldsc0(1,i)
 836 !        if (i.eq.10) then
 837 !          dsc_inv(i)=0.0D0
 838 !        else
 839 !          dsc_inv(i)=1.0D0/dsc(i)
 840 !        endif
 841 !      enddo
 842
 843       if (lprint) then
 844         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
 845         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
 846          'inertia','Pstok'
 847         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
 848         do i=1,ntyp
 849           write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
 850             vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
 851           do j=2,nbondterm(i)
 852             write (iout,'(13x,3f10.5)') &
 853               vbldsc0(j,i),aksc(j,i),abond0(j,i)
 854           enddo
 855         enddo
 856       endif
 857             do i=1,ntyp_molec(5)
 858              read(iion,*) msc(i,5),restok(i,5)
 859              print *,msc(i,5),restok(i,5)
 860             enddo
 861             ip(5)=0.2
 862 !            isc(5)=0.2
 863             read (iion,*) ncatprotparm
 864             allocate(catprm(ncatprotparm,4))
 865             do k=1,4
 866             read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
 867             enddo
 868             print *, catprm
 869 !            read (iion,*) (vcatprm(k),k=1,ncatprotpram)
 870 !----------------------------------------------------
 871       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
 872       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))      !(-ntyp:ntyp)
 873       allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
 874       allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
 875       allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
 876       allocate(gthet(3,-ntyp:ntyp))     !(3,-ntyp:ntyp)
 877
 878       a0thet(:)=0.0D0
 879       athet(:,:,:,:)=0.0D0
 880       bthet(:,:,:,:)=0.0D0
 881       polthet(:,:)=0.0D0
 882       gthet(:,:)=0.0D0
 883       theta0(:)=0.0D0
 884       sig0(:)=0.0D0
 885       sigc0(:)=0.0D0
 886       allocate(liptranene(ntyp))
 887 !C reading lipid parameters
 888       write (iout,*) "iliptranpar",iliptranpar
 889       call flush(iout)
 890        read(iliptranpar,*) pepliptran
 891        print *,pepliptran
 892        do i=1,ntyp
 893        read(iliptranpar,*) liptranene(i)
 894         print *,liptranene(i)
 895        enddo
 896        close(iliptranpar)
 897
 898 #ifdef CRYST_THETA
 899 !
 900 ! Read the parameters of the probability distribution/energy expression
 901 ! of the virtual-bond valence angles theta
 902 !
 903       do i=1,ntyp
 904         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
 905           (bthet(j,i,1,1),j=1,2)
 906         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 907         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 908         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 909         sigc0(i)=sigc0(i)**2
 910       enddo
 911       do i=1,ntyp
 912       athet(1,i,1,-1)=athet(1,i,1,1)
 913       athet(2,i,1,-1)=athet(2,i,1,1)
 914       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 915       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 916       athet(1,i,-1,1)=-athet(1,i,1,1)
 917       athet(2,i,-1,1)=-athet(2,i,1,1)
 918       bthet(1,i,-1,1)=bthet(1,i,1,1)
 919       bthet(2,i,-1,1)=bthet(2,i,1,1)
 920       enddo
 921       do i=-ntyp,-1
 922       a0thet(i)=a0thet(-i)
 923       athet(1,i,-1,-1)=athet(1,-i,1,1)
 924       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 925       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 926       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 927       athet(1,i,-1,1)=athet(1,-i,1,1)
 928       athet(2,i,-1,1)=-athet(2,-i,1,1)
 929       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 930       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 931       athet(1,i,1,-1)=-athet(1,-i,1,1)
 932       athet(2,i,1,-1)=athet(2,-i,1,1)
 933       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 934       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 935       theta0(i)=theta0(-i)
 936       sig0(i)=sig0(-i)
 937       sigc0(i)=sigc0(-i)
 938        do j=0,3
 939         polthet(j,i)=polthet(j,-i)
 940        enddo
 941        do j=1,3
 942          gthet(j,i)=gthet(j,-i)
 943        enddo
 944       enddo
 945
 946       close (ithep)
 947       if (lprint) then
 948       if (.not.LaTeX) then
 949         write (iout,'(a)') &
 950           'Parameters of the virtual-bond valence angles:'
 951         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
 952        '    ATHETA0   ','         A1   ','        A2    ',&
 953        '        B1    ','         B2   '
 954         do i=1,ntyp
 955           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
 956               a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
 957         enddo
 958         write (iout,'(/a/9x,5a/79(1h-))') &
 959        'Parameters of the expression for sigma(theta_c):',&
 960        '     ALPH0    ','      ALPH1   ','     ALPH2    ',&
 961        '     ALPH3    ','    SIGMA0C   '
 962         do i=1,ntyp
 963           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
 964             (polthet(j,i),j=0,3),sigc0(i)
 965         enddo
 966         write (iout,'(/a/9x,5a/79(1h-))') &
 967        'Parameters of the second gaussian:',&
 968        '    THETA0    ','     SIGMA0   ','        G1    ',&
 969        '        G2    ','         G3   '
 970         do i=1,ntyp
 971           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
 972              sig0(i),(gthet(j,i),j=1,3)
 973         enddo
 974        else
 975         write (iout,'(a)') &
 976           'Parameters of the virtual-bond valence angles:'
 977         write (iout,'(/a/9x,5a/79(1h-))') &
 978        'Coefficients of expansion',&
 979        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
 980        '   b1*10^1    ','    b2*10^1   '
 981         do i=1,ntyp
 982           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
 983               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
 984               (10*bthet(j,i,1,1),j=1,2)
 985         enddo
 986         write (iout,'(/a/9x,5a/79(1h-))') &
 987        'Parameters of the expression for sigma(theta_c):',&
 988        ' alpha0       ','  alph1       ',' alph2        ',&
 989        ' alhp3        ','   sigma0c    '
 990         do i=1,ntyp
 991           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
 992             (polthet(j,i),j=0,3),sigc0(i)
 993         enddo
 994         write (iout,'(/a/9x,5a/79(1h-))') &
 995        'Parameters of the second gaussian:',&
 996        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
 997        '        G2    ','   G3*10^1    '
 998         do i=1,ntyp
 999           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1000              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1001         enddo
1002       endif
1003       endif
1004 #else
1005 !
1006 ! Read the parameters of Utheta determined from ab initio surfaces
1007 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1008 !
1009       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1010         ntheterm3,nsingle,ndouble
1011       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1012
1013 !----------------------------------------------------
1014       allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1015       allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1016         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1017 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1018       allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1019         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1020 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1021 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1022       allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1023         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1024       allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1025         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1026       allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1027         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1028       allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1029         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1030 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1031 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1032       allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1033         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034       allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1035         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1036 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1037 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1038
1039       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1040       do i=-ntyp1,-1
1041         ithetyp(i)=-ithetyp(-i)
1042       enddo
1043
1044       aa0thet(:,:,:,:)=0.0d0
1045       aathet(:,:,:,:,:)=0.0d0
1046       bbthet(:,:,:,:,:,:)=0.0d0
1047       ccthet(:,:,:,:,:,:)=0.0d0
1048       ddthet(:,:,:,:,:,:)=0.0d0
1049       eethet(:,:,:,:,:,:)=0.0d0
1050       ffthet(:,:,:,:,:,:,:)=0.0d0
1051       ggthet(:,:,:,:,:,:,:)=0.0d0
1052
1053 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1054       do iblock=1,2
1055 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1056 ! VAR:1=non-glicyne non-proline 2=proline
1057 ! VAR:negative values for D-aminoacid
1058       do i=0,nthetyp
1059         do j=-nthetyp,nthetyp
1060           do k=-nthetyp,nthetyp
1061             read (ithep,'(6a)',end=111,err=111) res1
1062             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1063 ! VAR: aa0thet is variable describing the average value of Foureir
1064 ! VAR: expansion series
1065 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1066 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1067 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1068             read (ithep,*,end=111,err=111) &
1069               (aathet(l,i,j,k,iblock),l=1,ntheterm)
1070             read (ithep,*,end=111,err=111) &
1071              ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1072               (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1073               (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1074               (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1075               ll=1,ntheterm2)
1076             read (ithep,*,end=111,err=111) &
1077             (((ffthet(llll,lll,ll,i,j,k,iblock),&
1078                ffthet(lll,llll,ll,i,j,k,iblock),&
1079                ggthet(llll,lll,ll,i,j,k,iblock),&
1080                ggthet(lll,llll,ll,i,j,k,iblock),&
1081                llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1082           enddo
1083         enddo
1084       enddo
1085 !
1086 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1087 ! coefficients of theta-and-gamma-dependent terms are zero.
1088 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1089 ! RECOMENTDED AFTER VERSION 3.3)
1090 !      do i=1,nthetyp
1091 !        do j=1,nthetyp
1092 !          do l=1,ntheterm
1093 !            aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1094 !            aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1095 !          enddo
1096 !          aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1097 !          aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1098 !        enddo
1099 !        do l=1,ntheterm
1100 !          aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1101 !        enddo
1102 !        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1103 !      enddo
1104 !      enddo
1105 ! AND COMMENT THE LOOPS BELOW
1106       do i=1,nthetyp
1107         do j=1,nthetyp
1108           do l=1,ntheterm
1109             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1110             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1111           enddo
1112           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1113           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1114         enddo
1115         do l=1,ntheterm
1116           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1117         enddo
1118         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1119       enddo
1120       enddo       !iblock
1121
1122 ! TILL HERE
1123 ! Substitution for D aminoacids from symmetry.
1124       do iblock=1,2
1125       do i=-nthetyp,0
1126         do j=-nthetyp,nthetyp
1127           do k=-nthetyp,nthetyp
1128            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1129            do l=1,ntheterm
1130            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1131            enddo
1132            do ll=1,ntheterm2
1133             do lll=1,nsingle
1134             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1135             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1136             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1137             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1138             enddo
1139           enddo
1140           do ll=1,ntheterm3
1141            do lll=2,ndouble
1142             do llll=1,lll-1
1143             ffthet(llll,lll,ll,i,j,k,iblock)= &
1144             ffthet(llll,lll,ll,-i,-j,-k,iblock)
1145             ffthet(lll,llll,ll,i,j,k,iblock)= &
1146             ffthet(lll,llll,ll,-i,-j,-k,iblock)
1147             ggthet(llll,lll,ll,i,j,k,iblock)= &
1148             -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1149             ggthet(lll,llll,ll,i,j,k,iblock)= &
1150             -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1151             enddo !ll
1152            enddo  !lll
1153           enddo   !llll
1154          enddo    !k
1155         enddo     !j
1156        enddo      !i
1157       enddo       !iblock
1158 !
1159 ! Control printout of the coefficients of virtual-bond-angle potentials
1160 !
1161       if (lprint) then
1162         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1163         do iblock=1,2
1164         do i=1,nthetyp+1
1165           do j=1,nthetyp+1
1166             do k=1,nthetyp+1
1167               write (iout,'(//4a)') &
1168                'Type ',onelett(i),onelett(j),onelett(k)
1169               write (iout,'(//a,10x,a)') " l","a[l]"
1170               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1171               write (iout,'(i2,1pe15.5)') &
1172                  (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1173             do l=1,ntheterm2
1174               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1175                 "b",l,"c",l,"d",l,"e",l
1176               do m=1,nsingle
1177                 write (iout,'(i2,4(1pe15.5))') m,&
1178                 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1179                 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1180               enddo
1181             enddo
1182             do l=1,ntheterm3
1183               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1184                 "f+",l,"f-",l,"g+",l,"g-",l
1185               do m=2,ndouble
1186                 do n=1,m-1
1187                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1188                     ffthet(n,m,l,i,j,k,iblock),&
1189                     ffthet(m,n,l,i,j,k,iblock),&
1190                     ggthet(n,m,l,i,j,k,iblock),&
1191                     ggthet(m,n,l,i,j,k,iblock)
1192                 enddo   !n
1193               enddo     !m
1194             enddo       !l
1195           enddo         !k
1196         enddo           !j
1197       enddo             !i
1198       enddo
1199       call flush(iout)
1200       endif
1201       write (2,*) "Start reading THETA_PDB",ithep_pdb
1202       do i=1,ntyp
1203 !      write (2,*) 'i=',i
1204         read (ithep_pdb,*,err=111,end=111) &
1205            a0thet(i),(athet(j,i,1,1),j=1,2),&
1206           (bthet(j,i,1,1),j=1,2)
1207         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1208         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1209         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1210         sigc0(i)=sigc0(i)**2
1211       enddo
1212       do i=1,ntyp
1213       athet(1,i,1,-1)=athet(1,i,1,1)
1214       athet(2,i,1,-1)=athet(2,i,1,1)
1215       bthet(1,i,1,-1)=-bthet(1,i,1,1)
1216       bthet(2,i,1,-1)=-bthet(2,i,1,1)
1217       athet(1,i,-1,1)=-athet(1,i,1,1)
1218       athet(2,i,-1,1)=-athet(2,i,1,1)
1219       bthet(1,i,-1,1)=bthet(1,i,1,1)
1220       bthet(2,i,-1,1)=bthet(2,i,1,1)
1221       enddo
1222       do i=-ntyp,-1
1223       a0thet(i)=a0thet(-i)
1224       athet(1,i,-1,-1)=athet(1,-i,1,1)
1225       athet(2,i,-1,-1)=-athet(2,-i,1,1)
1226       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1227       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1228       athet(1,i,-1,1)=athet(1,-i,1,1)
1229       athet(2,i,-1,1)=-athet(2,-i,1,1)
1230       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1231       bthet(2,i,-1,1)=bthet(2,-i,1,1)
1232       athet(1,i,1,-1)=-athet(1,-i,1,1)
1233       athet(2,i,1,-1)=athet(2,-i,1,1)
1234       bthet(1,i,1,-1)=bthet(1,-i,1,1)
1235       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1236       theta0(i)=theta0(-i)
1237       sig0(i)=sig0(-i)
1238       sigc0(i)=sigc0(-i)
1239        do j=0,3
1240         polthet(j,i)=polthet(j,-i)
1241        enddo
1242        do j=1,3
1243          gthet(j,i)=gthet(j,-i)
1244        enddo
1245       enddo
1246       write (2,*) "End reading THETA_PDB"
1247       close (ithep_pdb)
1248 #endif
1249       close(ithep)
1250 !--------------- Reading theta parameters for nucleic acid-------
1251       read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1252       ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1253       nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1254       allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1255       allocate(aa0thet_nucl(nthetyp_nucl+1,&
1256         nthetyp_nucl+1,nthetyp_nucl+1))
1257 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1258       allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1259         nthetyp_nucl+1,nthetyp_nucl+1))
1260 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1261 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1262       allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1263         nthetyp_nucl+1,nthetyp_nucl+1))
1264       allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1265         nthetyp_nucl+1,nthetyp_nucl+1))
1266       allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1267         nthetyp_nucl+1,nthetyp_nucl+1))
1268       allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1269         nthetyp_nucl+1,nthetyp_nucl+1))
1270 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1271 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1272       allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1273          nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1274       allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1275          nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1276
1277 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1278 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1279
1280       read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1281
1282       aa0thet_nucl(:,:,:)=0.0d0
1283       aathet_nucl(:,:,:,:)=0.0d0
1284       bbthet_nucl(:,:,:,:,:)=0.0d0
1285       ccthet_nucl(:,:,:,:,:)=0.0d0
1286       ddthet_nucl(:,:,:,:,:)=0.0d0
1287       eethet_nucl(:,:,:,:,:)=0.0d0
1288       ffthet_nucl(:,:,:,:,:,:)=0.0d0
1289       ggthet_nucl(:,:,:,:,:,:)=0.0d0
1290
1291       do i=1,nthetyp_nucl
1292         do j=1,nthetyp_nucl
1293           do k=1,nthetyp_nucl
1294             read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1295             read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1296             read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1297             read (ithep_nucl,*,end=111,err=111) &
1298             (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1299             (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1300             (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1301             (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1302             read (ithep_nucl,*,end=111,err=111) &
1303            (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1304               ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1305               llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1306           enddo
1307         enddo
1308       enddo
1309
1310 !-------------------------------------------
1311       allocate(nlob(ntyp1)) !(ntyp1)
1312       allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1313       allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1314       allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1315
1316       bsc(:,:)=0.0D0
1317       nlob(:)=0
1318       nlob(:)=0
1319       dsc(:)=0.0D0
1320       censc(:,:,:)=0.0D0
1321       gaussc(:,:,:,:)=0.0D0
1322
1323 #ifdef CRYST_SC
1324 !
1325 ! Read the parameters of the probability distribution/energy expression
1326 ! of the side chains.
1327 !
1328       do i=1,ntyp
1329         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1330         if (i.eq.10) then
1331           dsc_inv(i)=0.0D0
1332         else
1333           dsc_inv(i)=1.0D0/dsc(i)
1334         endif
1335         if (i.ne.10) then
1336         do j=1,nlob(i)
1337           do k=1,3
1338             do l=1,3
1339               blower(l,k,j)=0.0D0
1340             enddo
1341           enddo
1342         enddo
1343         bsc(1,i)=0.0D0
1344         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1345           ((blower(k,l,1),l=1,k),k=1,3)
1346         censc(1,1,-i)=censc(1,1,i)
1347         censc(2,1,-i)=censc(2,1,i)
1348         censc(3,1,-i)=-censc(3,1,i)
1349         do j=2,nlob(i)
1350           read (irotam,*,end=112,err=112) bsc(j,i)
1351           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1352                                        ((blower(k,l,j),l=1,k),k=1,3)
1353         censc(1,j,-i)=censc(1,j,i)
1354         censc(2,j,-i)=censc(2,j,i)
1355         censc(3,j,-i)=-censc(3,j,i)
1356 ! BSC is amplitude of Gaussian
1357         enddo
1358         do j=1,nlob(i)
1359           do k=1,3
1360             do l=1,k
1361               akl=0.0D0
1362               do m=1,3
1363                 akl=akl+blower(k,m,j)*blower(l,m,j)
1364               enddo
1365               gaussc(k,l,j,i)=akl
1366               gaussc(l,k,j,i)=akl
1367              if (((k.eq.3).and.(l.ne.3)) &
1368               .or.((l.eq.3).and.(k.ne.3))) then
1369                 gaussc(k,l,j,-i)=-akl
1370                 gaussc(l,k,j,-i)=-akl
1371               else
1372                 gaussc(k,l,j,-i)=akl
1373                 gaussc(l,k,j,-i)=akl
1374               endif
1375             enddo
1376           enddo
1377         enddo
1378         endif
1379       enddo
1380       close (irotam)
1381       if (lprint) then
1382         write (iout,'(/a)') 'Parameters of side-chain local geometry'
1383         do i=1,ntyp
1384           nlobi=nlob(i)
1385           if (nlobi.gt.0) then
1386             if (LaTeX) then
1387               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1388                ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1389                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1390                                    'log h',(bsc(j,i),j=1,nlobi)
1391                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1392               'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1393               do k=1,3
1394                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1395                        ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1396               enddo
1397             else
1398               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1399               write (iout,'(a,f10.4,4(16x,f10.4))') &
1400                                    'Center  ',(bsc(j,i),j=1,nlobi)
1401               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1402                  j=1,nlobi)
1403               write (iout,'(a)')
1404             endif
1405           endif
1406         enddo
1407       endif
1408 #else
1409 !
1410 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1411 ! added by Urszula Kozlowska 07/11/2007
1412 !
1413 !el Maximum number of SC local term fitting function coefficiants
1414 !el      integer,parameter :: maxsccoef=65
1415
1416       allocate(sc_parmin(65,ntyp))      !(maxsccoef,ntyp)
1417
1418       do i=1,ntyp
1419         read (irotam,*,end=112,err=112)
1420        if (i.eq.10) then
1421          read (irotam,*,end=112,err=112)
1422        else
1423          do j=1,65
1424            read(irotam,*,end=112,err=112) sc_parmin(j,i)
1425          enddo
1426        endif
1427       enddo
1428 !---------reading nucleic acid parameters for rotamers-------------------
1429       allocate(sc_parmin_nucl(9,ntyp_molec(2)))      !(maxsccoef,ntyp)
1430       do i=1,ntyp_molec(2)
1431         read (irotam_nucl,*,end=112,err=112)
1432         do j=1,9
1433           read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1434         enddo
1435       enddo
1436       close(irotam_nucl)
1437       if (lprint) then
1438         write (iout,*)
1439         write (iout,*) "Base rotamer parameters"
1440         do i=1,ntyp_molec(2)
1441           write (iout,'(a)') restyp(i,2)
1442           write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1443         enddo
1444       endif
1445
1446 !
1447 ! Read the parameters of the probability distribution/energy expression
1448 ! of the side chains.
1449 !
1450       write (2,*) "Start reading ROTAM_PDB"
1451       do i=1,ntyp
1452         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1453         if (i.eq.10) then
1454           dsc_inv(i)=0.0D0
1455         else
1456           dsc_inv(i)=1.0D0/dsc(i)
1457         endif
1458         if (i.ne.10) then
1459         do j=1,nlob(i)
1460           do k=1,3
1461             do l=1,3
1462               blower(l,k,j)=0.0D0
1463             enddo
1464           enddo
1465         enddo
1466         bsc(1,i)=0.0D0
1467         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1468           ((blower(k,l,1),l=1,k),k=1,3)
1469         do j=2,nlob(i)
1470           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1471           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1472                                        ((blower(k,l,j),l=1,k),k=1,3)
1473         enddo
1474         do j=1,nlob(i)
1475           do k=1,3
1476             do l=1,k
1477               akl=0.0D0
1478               do m=1,3
1479                 akl=akl+blower(k,m,j)*blower(l,m,j)
1480               enddo
1481               gaussc(k,l,j,i)=akl
1482               gaussc(l,k,j,i)=akl
1483             enddo
1484           enddo
1485         enddo
1486         endif
1487       enddo
1488       close (irotam_pdb)
1489       write (2,*) "End reading ROTAM_PDB"
1490 #endif
1491       close(irotam)
1492
1493 #ifdef CRYST_TOR
1494 !
1495 ! Read torsional parameters in old format
1496 !
1497       allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1498
1499       read (itorp,*,end=113,err=113) ntortyp,nterm_old
1500       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1501       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1502
1503 !el from energy module--------
1504       allocate(v1(nterm_old,ntortyp,ntortyp))
1505       allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1506 !el---------------------------
1507       do i=1,ntortyp
1508         do j=1,ntortyp
1509           read (itorp,'(a)')
1510           do k=1,nterm_old
1511             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1512           enddo
1513         enddo
1514       enddo
1515       close (itorp)
1516       if (lprint) then
1517         write (iout,'(/a/)') 'Torsional constants:'
1518         do i=1,ntortyp
1519           do j=1,ntortyp
1520             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1521             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1522           enddo
1523         enddo
1524       endif
1525 #else
1526 !
1527 ! Read torsional parameters
1528 !
1529       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1530       read (itorp,*,end=113,err=113) ntortyp
1531 !el from energy module---------
1532       allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1533       allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1534
1535       allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1536       allocate(vlor2(maxlor,ntortyp,ntortyp))
1537       allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1538       allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1539
1540       allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1541       allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1542 !el---------------------------
1543       nterm(:,:,:)=0
1544       nlor(:,:,:)=0
1545 !el---------------------------
1546
1547       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1548       do i=-ntyp,-1
1549        itortyp(i)=-itortyp(-i)
1550       enddo
1551       itortyp(ntyp1)=ntortyp
1552       itortyp(-ntyp1)=-ntortyp
1553       do iblock=1,2
1554       write (iout,*) 'ntortyp',ntortyp
1555       do i=0,ntortyp-1
1556         do j=-ntortyp+1,ntortyp-1
1557           read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1558                 nlor(i,j,iblock)
1559           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1560           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1561           v0ij=0.0d0
1562           si=-1.0d0
1563           do k=1,nterm(i,j,iblock)
1564             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1565             v2(k,i,j,iblock)
1566             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1567             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1568             v0ij=v0ij+si*v1(k,i,j,iblock)
1569             si=-si
1570 !         write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1571 !         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1572 !      v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1573           enddo
1574           do k=1,nlor(i,j,iblock)
1575             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1576               vlor2(k,i,j),vlor3(k,i,j)
1577             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1578           enddo
1579           v0(i,j,iblock)=v0ij
1580           v0(-i,-j,iblock)=v0ij
1581         enddo
1582       enddo
1583       enddo
1584       close (itorp)
1585       if (lprint) then
1586         write (iout,'(/a/)') 'Torsional constants:'
1587         do iblock=1,2
1588         do i=-ntortyp,ntortyp
1589           do j=-ntortyp,ntortyp
1590             write (iout,*) 'ityp',i,' jtyp',j
1591             write (iout,*) 'Fourier constants'
1592             do k=1,nterm(i,j,iblock)
1593               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1594               v2(k,i,j,iblock)
1595             enddo
1596             write (iout,*) 'Lorenz constants'
1597             do k=1,nlor(i,j,iblock)
1598               write (iout,'(3(1pe15.5))') &
1599                vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1600             enddo
1601           enddo
1602         enddo
1603         enddo
1604       endif
1605 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1606 !
1607 ! 6/23/01 Read parameters for double torsionals
1608 !
1609 !el from energy module------------
1610       allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1611       allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1612 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1613       allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1614       allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1615         !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1616       allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1617       allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1618         !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1619 !---------------------------------
1620
1621       do iblock=1,2
1622       do i=0,ntortyp-1
1623         do j=-ntortyp+1,ntortyp-1
1624           do k=-ntortyp+1,ntortyp-1
1625             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1626 !              write (iout,*) "OK onelett",
1627 !     &         i,j,k,t1,t2,t3
1628
1629             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1630               .or. t3.ne.toronelet(k)) then
1631              write (iout,*) "Error in double torsional parameter file",&
1632                i,j,k,t1,t2,t3
1633 #ifdef MPI
1634               call MPI_Finalize(Ierror)
1635 #endif
1636                stop "Error in double torsional parameter file"
1637             endif
1638            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1639                ntermd_2(i,j,k,iblock)
1640             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1641             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1642             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1643                ntermd_1(i,j,k,iblock))
1644             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1645                ntermd_1(i,j,k,iblock))
1646             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1647                ntermd_1(i,j,k,iblock))
1648             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1649                ntermd_1(i,j,k,iblock))
1650 ! Martix of D parameters for one dimesional foureir series
1651             do l=1,ntermd_1(i,j,k,iblock)
1652              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1653              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1654              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1655              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1656 !            write(iout,*) "whcodze" ,
1657 !     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1658             enddo
1659             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1660                v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1661                v2s(m,l,i,j,k,iblock),&
1662                m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1663 ! Martix of D parameters for two dimesional fourier series
1664             do l=1,ntermd_2(i,j,k,iblock)
1665              do m=1,l-1
1666              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1667              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1668              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1669              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1670              enddo!m
1671             enddo!l
1672           enddo!k
1673         enddo!j
1674       enddo!i
1675       enddo!iblock
1676       if (lprint) then
1677       write (iout,*)
1678       write (iout,*) 'Constants for double torsionals'
1679       do iblock=1,2
1680       do i=0,ntortyp-1
1681         do j=-ntortyp+1,ntortyp-1
1682           do k=-ntortyp+1,ntortyp-1
1683             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1684               ' nsingle',ntermd_1(i,j,k,iblock),&
1685               ' ndouble',ntermd_2(i,j,k,iblock)
1686             write (iout,*)
1687             write (iout,*) 'Single angles:'
1688             do l=1,ntermd_1(i,j,k,iblock)
1689               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1690                  v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1691                  v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1692                  v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1693             enddo
1694             write (iout,*)
1695             write (iout,*) 'Pairs of angles:'
1696             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1697             do l=1,ntermd_2(i,j,k,iblock)
1698               write (iout,'(i5,20f10.5)') &
1699                l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1700             enddo
1701             write (iout,*)
1702             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1703             do l=1,ntermd_2(i,j,k,iblock)
1704               write (iout,'(i5,20f10.5)') &
1705                l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1706                (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1707             enddo
1708             write (iout,*)
1709           enddo
1710         enddo
1711       enddo
1712       enddo
1713       endif
1714 #endif
1715       allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1716       read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1717 !      print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1718 !el from energy module---------
1719       allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1720       allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1721
1722       allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1723       allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1724       allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1725       allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1726
1727       allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1728       allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1729 !el---------------------------
1730       nterm_nucl(:,:)=0
1731       nlor_nucl(:,:)=0
1732 !el--------------------
1733       read (itorp_nucl,*,end=113,err=113) &
1734         (itortyp_nucl(i),i=1,ntyp_molec(2))
1735 !        print *,itortyp_nucl(:)
1736 !c      write (iout,*) 'ntortyp',ntortyp
1737       do i=1,ntortyp_nucl
1738         do j=1,ntortyp_nucl
1739           read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1740 !           print *,nterm_nucl(i,j),nlor_nucl(i,j)
1741           v0ij=0.0d0
1742           si=-1.0d0
1743           do k=1,nterm_nucl(i,j)
1744             read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1745             v0ij=v0ij+si*v1_nucl(k,i,j)
1746             si=-si
1747           enddo
1748           do k=1,nlor_nucl(i,j)
1749             read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1750               vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1751             v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1752           enddo
1753           v0_nucl(i,j)=v0ij
1754         enddo
1755       enddo
1756
1757 ! Read of Side-chain backbone correlation parameters
1758 ! Modified 11 May 2012 by Adasko
1759 !CC
1760 !
1761       read (isccor,*,end=119,err=119) nsccortyp
1762
1763 !el from module energy-------------
1764       allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1765       allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1766       allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1767       allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp))   !(maxterm_sccor,20,20)
1768 !-----------------------------------
1769 #ifdef SCCORPDB
1770 !el from module energy-------------
1771       allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1772
1773       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1774       do i=-ntyp,-1
1775         isccortyp(i)=-isccortyp(-i)
1776       enddo
1777       iscprol=isccortyp(20)
1778 !      write (iout,*) 'ntortyp',ntortyp
1779       maxinter=3
1780 !c maxinter is maximum interaction sites
1781 !el from module energy---------
1782       allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1783       allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1784                -nsccortyp:nsccortyp))
1785       allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1786                -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1787       allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1788                -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1789 !-----------------------------------
1790       do l=1,maxinter
1791       do i=1,nsccortyp
1792         do j=1,nsccortyp
1793           read (isccor,*,end=119,err=119) &
1794       nterm_sccor(i,j),nlor_sccor(i,j)
1795           v0ijsccor=0.0d0
1796           v0ijsccor1=0.0d0
1797           v0ijsccor2=0.0d0
1798           v0ijsccor3=0.0d0
1799           si=-1.0d0
1800           nterm_sccor(-i,j)=nterm_sccor(i,j)
1801           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1802           nterm_sccor(i,-j)=nterm_sccor(i,j)
1803           do k=1,nterm_sccor(i,j)
1804             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1805            v2sccor(k,l,i,j)
1806             if (j.eq.iscprol) then
1807              if (i.eq.isccortyp(10)) then
1808              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1809              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1810              else
1811              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1812                               +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1813              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1814                               +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1815              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1816              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1817              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1818              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1819              endif
1820             else
1821              if (i.eq.isccortyp(10)) then
1822              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1823              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1824              else
1825                if (j.eq.isccortyp(10)) then
1826              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1827              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1828                else
1829              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1830              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1831              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1832              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1833              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1834              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1835                 endif
1836                endif
1837             endif
1838             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1839             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1840             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1841             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1842             si=-si
1843           enddo
1844           do k=1,nlor_sccor(i,j)
1845             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1846               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1847             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1848       (1+vlor3sccor(k,i,j)**2)
1849           enddo
1850           v0sccor(l,i,j)=v0ijsccor
1851           v0sccor(l,-i,j)=v0ijsccor1
1852           v0sccor(l,i,-j)=v0ijsccor2
1853           v0sccor(l,-i,-j)=v0ijsccor3
1854         enddo
1855       enddo
1856       enddo
1857       close (isccor)
1858 #else
1859 !el from module energy-------------
1860       allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1861
1862       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1863 !      write (iout,*) 'ntortyp',ntortyp
1864       maxinter=3
1865 !c maxinter is maximum interaction sites
1866 !el from module energy---------
1867       allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1868       allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1869       allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1870       allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1871 !-----------------------------------
1872       do l=1,maxinter
1873       do i=1,nsccortyp
1874         do j=1,nsccortyp
1875           read (isccor,*,end=119,err=119) &
1876        nterm_sccor(i,j),nlor_sccor(i,j)
1877           v0ijsccor=0.0d0
1878           si=-1.0d0
1879
1880           do k=1,nterm_sccor(i,j)
1881             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1882            v2sccor(k,l,i,j)
1883             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1884             si=-si
1885           enddo
1886           do k=1,nlor_sccor(i,j)
1887             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1888               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1889             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1890       (1+vlor3sccor(k,i,j)**2)
1891           enddo
1892           v0sccor(l,i,j)=v0ijsccor !el ,iblock
1893         enddo
1894       enddo
1895       enddo
1896       close (isccor)
1897
1898 #endif
1899       if (lprint) then
1900         write (iout,'(/a/)') 'Torsional constants:'
1901         do i=1,nsccortyp
1902           do j=1,nsccortyp
1903             write (iout,*) 'ityp',i,' jtyp',j
1904             write (iout,*) 'Fourier constants'
1905             do k=1,nterm_sccor(i,j)
1906       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1907             enddo
1908             write (iout,*) 'Lorenz constants'
1909             do k=1,nlor_sccor(i,j)
1910               write (iout,'(3(1pe15.5))') &
1911                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1912             enddo
1913           enddo
1914         enddo
1915       endif
1916
1917 !
1918 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1919 !         interaction energy of the Gly, Ala, and Pro prototypes.
1920 !
1921
1922       if (lprint) then
1923         write (iout,*)
1924         write (iout,*) "Coefficients of the cumulants"
1925       endif
1926       read (ifourier,*) nloctyp
1927 !write(iout,*) "nloctyp",nloctyp
1928 !el from module energy-------
1929       allocate(b1(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1930       allocate(b2(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1931       allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1932       allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1933       allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1934       allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1935       allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1936       allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1937 ! el
1938         b1(1,:)=0.0d0
1939         b1(2,:)=0.0d0
1940 !--------------------------------
1941
1942       do i=0,nloctyp-1
1943         read (ifourier,*,end=115,err=115)
1944         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1945         if (lprint) then
1946           write (iout,*) 'Type',i
1947           write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1948         endif
1949         B1(1,i)  = b(3)
1950         B1(2,i)  = b(5)
1951         B1(1,-i) = b(3)
1952         B1(2,-i) = -b(5)
1953 !        b1(1,i)=0.0d0
1954 !        b1(2,i)=0.0d0
1955         B1tilde(1,i) = b(3)
1956         B1tilde(2,i) =-b(5)
1957         B1tilde(1,-i) =-b(3)
1958         B1tilde(2,-i) =b(5)
1959 !        b1tilde(1,i)=0.0d0
1960 !        b1tilde(2,i)=0.0d0
1961         B2(1,i)  = b(2)
1962         B2(2,i)  = b(4)
1963         B2(1,-i)  =b(2)
1964         B2(2,-i)  =-b(4)
1965
1966 !        b2(1,i)=0.0d0
1967 !        b2(2,i)=0.0d0
1968         CC(1,1,i)= b(7)
1969         CC(2,2,i)=-b(7)
1970         CC(2,1,i)= b(9)
1971         CC(1,2,i)= b(9)
1972         CC(1,1,-i)= b(7)
1973         CC(2,2,-i)=-b(7)
1974         CC(2,1,-i)=-b(9)
1975         CC(1,2,-i)=-b(9)
1976 !        CC(1,1,i)=0.0d0
1977 !        CC(2,2,i)=0.0d0
1978 !        CC(2,1,i)=0.0d0
1979 !        CC(1,2,i)=0.0d0
1980         Ctilde(1,1,i)=b(7)
1981         Ctilde(1,2,i)=b(9)
1982         Ctilde(2,1,i)=-b(9)
1983         Ctilde(2,2,i)=b(7)
1984         Ctilde(1,1,-i)=b(7)
1985         Ctilde(1,2,-i)=-b(9)
1986         Ctilde(2,1,-i)=b(9)
1987         Ctilde(2,2,-i)=b(7)
1988
1989 !        Ctilde(1,1,i)=0.0d0
1990 !        Ctilde(1,2,i)=0.0d0
1991 !        Ctilde(2,1,i)=0.0d0
1992 !        Ctilde(2,2,i)=0.0d0
1993         DD(1,1,i)= b(6)
1994         DD(2,2,i)=-b(6)
1995         DD(2,1,i)= b(8)
1996         DD(1,2,i)= b(8)
1997         DD(1,1,-i)= b(6)
1998         DD(2,2,-i)=-b(6)
1999         DD(2,1,-i)=-b(8)
2000         DD(1,2,-i)=-b(8)
2001 !        DD(1,1,i)=0.0d0
2002 !        DD(2,2,i)=0.0d0
2003 !        DD(2,1,i)=0.0d0
2004 !        DD(1,2,i)=0.0d0
2005         Dtilde(1,1,i)=b(6)
2006         Dtilde(1,2,i)=b(8)
2007         Dtilde(2,1,i)=-b(8)
2008         Dtilde(2,2,i)=b(6)
2009         Dtilde(1,1,-i)=b(6)
2010         Dtilde(1,2,-i)=-b(8)
2011         Dtilde(2,1,-i)=b(8)
2012         Dtilde(2,2,-i)=b(6)
2013
2014 !        Dtilde(1,1,i)=0.0d0
2015 !        Dtilde(1,2,i)=0.0d0
2016 !        Dtilde(2,1,i)=0.0d0
2017 !        Dtilde(2,2,i)=0.0d0
2018         EE(1,1,i)= b(10)+b(11)
2019         EE(2,2,i)=-b(10)+b(11)
2020         EE(2,1,i)= b(12)-b(13)
2021         EE(1,2,i)= b(12)+b(13)
2022         EE(1,1,-i)= b(10)+b(11)
2023         EE(2,2,-i)=-b(10)+b(11)
2024         EE(2,1,-i)=-b(12)+b(13)
2025         EE(1,2,-i)=-b(12)-b(13)
2026
2027 !        ee(1,1,i)=1.0d0
2028 !        ee(2,2,i)=1.0d0
2029 !        ee(2,1,i)=0.0d0
2030 !        ee(1,2,i)=0.0d0
2031 !        ee(2,1,i)=ee(1,2,i)
2032       enddo
2033       if (lprint) then
2034       do i=1,nloctyp
2035         write (iout,*) 'Type',i
2036         write (iout,*) 'B1'
2037         write(iout,*) B1(1,i),B1(2,i)
2038         write (iout,*) 'B2'
2039         write(iout,*) B2(1,i),B2(2,i)
2040         write (iout,*) 'CC'
2041         do j=1,2
2042           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2043         enddo
2044         write(iout,*) 'DD'
2045         do j=1,2
2046           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2047         enddo
2048         write(iout,*) 'EE'
2049         do j=1,2
2050           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2051         enddo
2052       enddo
2053       endif
2054 !
2055 ! Read electrostatic-interaction parameters
2056 !
2057
2058       if (lprint) then
2059         write (iout,*)
2060         write (iout,'(/a)') 'Electrostatic interaction constants:'
2061         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2062                   'IT','JT','APP','BPP','AEL6','AEL3'
2063       endif
2064       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2065       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2066       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2067       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2068       close (ielep)
2069       do i=1,2
2070         do j=1,2
2071         rri=rpp(i,j)**6
2072         app (i,j)=epp(i,j)*rri*rri
2073         bpp (i,j)=-2.0D0*epp(i,j)*rri
2074         ael6(i,j)=elpp6(i,j)*4.2D0**6
2075         ael3(i,j)=elpp3(i,j)*4.2D0**3
2076 !        lprint=.true.
2077         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2078                           ael6(i,j),ael3(i,j)
2079 !        lprint=.false.
2080         enddo
2081       enddo
2082 !
2083 ! Read side-chain interaction parameters.
2084 !
2085 !el from module energy - COMMON.INTERACT-------
2086       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2087       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2088       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2089
2090       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2091       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2092       allocate(epslip(ntyp,ntyp))
2093       augm(:,:)=0.0D0
2094       chip(:)=0.0D0
2095       alp(:)=0.0D0
2096       sigma0(:)=0.0D0
2097       sigii(:)=0.0D0
2098       rr0(:)=0.0D0
2099
2100 !--------------------------------
2101
2102       read (isidep,*,end=117,err=117) ipot,expon
2103       if (ipot.lt.1 .or. ipot.gt.5) then
2104         write (iout,'(2a)') 'Error while reading SC interaction',&
2105                      'potential file - unknown potential type.'
2106 #ifdef MPI
2107         call MPI_Finalize(Ierror)
2108 #endif
2109         stop
2110       endif
2111       expon2=expon/2
2112       if(me.eq.king) &
2113        write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2114        ', exponents are ',expon,2*expon
2115 !      goto (10,20,30,30,40) ipot
2116       select case(ipot)
2117 !----------------------- LJ potential ---------------------------------
2118        case (1)
2119 !   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2120          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2121            (sigma0(i),i=1,ntyp)
2122         if (lprint) then
2123           write (iout,'(/a/)') 'Parameters of the LJ potential:'
2124           write (iout,'(a/)') 'The epsilon array:'
2125           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2126           write (iout,'(/a)') 'One-body parameters:'
2127           write (iout,'(a,4x,a)') 'residue','sigma'
2128           write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2129         endif
2130 !      goto 50
2131 !----------------------- LJK potential --------------------------------
2132        case(2)
2133 !   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2134          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2135           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2136         if (lprint) then
2137           write (iout,'(/a/)') 'Parameters of the LJK potential:'
2138           write (iout,'(a/)') 'The epsilon array:'
2139           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2140           write (iout,'(/a)') 'One-body parameters:'
2141           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
2142           write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2143                 i=1,ntyp)
2144         endif
2145 !      goto 50
2146 !---------------------- GB or BP potential -----------------------------
2147        case(3:4)
2148 !   30 do i=1,ntyp
2149         do i=1,ntyp
2150          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2151         enddo
2152         read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2153         read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2154         read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2155         read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2156         do i=1,ntyp
2157          read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2158         enddo
2159
2160 ! For the GB potential convert sigma'**2 into chi'
2161         if (ipot.eq.4) then
2162           do i=1,ntyp
2163             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2164           enddo
2165         endif
2166         if (lprint) then
2167           write (iout,'(/a/)') 'Parameters of the BP potential:'
2168           write (iout,'(a/)') 'The epsilon array:'
2169           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2170           write (iout,'(/a)') 'One-body parameters:'
2171           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
2172                '    chip  ','    alph  '
2173           write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2174                              chip(i),alp(i),i=1,ntyp)
2175         endif
2176 !      goto 50
2177 !--------------------- GBV potential -----------------------------------
2178        case(5)
2179 !   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2180         read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2181           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2182           (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2183         if (lprint) then
2184           write (iout,'(/a/)') 'Parameters of the GBV potential:'
2185           write (iout,'(a/)') 'The epsilon array:'
2186           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2187           write (iout,'(/a)') 'One-body parameters:'
2188           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
2189               's||/s_|_^2','    chip  ','    alph  '
2190           write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2191                    sigii(i),chip(i),alp(i),i=1,ntyp)
2192         endif
2193        case default
2194         write(iout,*)"Wrong ipot"
2195 !   50 continue
2196       end select
2197       continue
2198       close (isidep)
2199
2200 !-----------------------------------------------------------------------
2201 ! Calculate the "working" parameters of SC interactions.
2202
2203 !el from module energy - COMMON.INTERACT-------
2204       allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2205       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2206       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2207       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2208         dcavtub(ntyp1))
2209       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2210         tubetranene(ntyp1))
2211       aa_aq(:,:)=0.0D0
2212       bb_aq(:,:)=0.0D0
2213       aa_lip(:,:)=0.0D0
2214       bb_lip(:,:)=0.0D0
2215       chi(:,:)=0.0D0
2216       sigma(:,:)=0.0D0
2217       r0(:,:)=0.0D0
2218       acavtub(:)=0.0d0
2219       bcavtub(:)=0.0d0
2220       ccavtub(:)=0.0d0
2221       dcavtub(:)=0.0d0
2222       sc_aa_tube_par(:)=0.0d0
2223       sc_bb_tube_par(:)=0.0d0
2224
2225 !--------------------------------
2226
2227       do i=2,ntyp
2228         do j=1,i-1
2229           eps(i,j)=eps(j,i)
2230           epslip(i,j)=epslip(j,i)
2231         enddo
2232       enddo
2233       do i=1,ntyp
2234         do j=i,ntyp
2235           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2236           sigma(j,i)=sigma(i,j)
2237           rs0(i,j)=dwa16*sigma(i,j)
2238           rs0(j,i)=rs0(i,j)
2239         enddo
2240       enddo
2241       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2242        'Working parameters of the SC interactions:',&
2243        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
2244        '  chi1   ','   chi2   '
2245       do i=1,ntyp
2246         do j=i,ntyp
2247           epsij=eps(i,j)
2248           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2249             rrij=sigma(i,j)
2250           else
2251             rrij=rr0(i)+rr0(j)
2252           endif
2253           r0(i,j)=rrij
2254           r0(j,i)=rrij
2255           rrij=rrij**expon
2256           epsij=eps(i,j)
2257           sigeps=dsign(1.0D0,epsij)
2258           epsij=dabs(epsij)
2259           aa_aq(i,j)=epsij*rrij*rrij
2260           bb_aq(i,j)=-sigeps*epsij*rrij
2261           aa_aq(j,i)=aa_aq(i,j)
2262           bb_aq(j,i)=bb_aq(i,j)
2263           epsijlip=epslip(i,j)
2264           sigeps=dsign(1.0D0,epsijlip)
2265           epsijlip=dabs(epsijlip)
2266           aa_lip(i,j)=epsijlip*rrij*rrij
2267           bb_lip(i,j)=-sigeps*epsijlip*rrij
2268           aa_lip(j,i)=aa_lip(i,j)
2269           bb_lip(j,i)=bb_lip(i,j)
2270 !C          write(iout,*) aa_lip
2271           if (ipot.gt.2) then
2272             sigt1sq=sigma0(i)**2
2273             sigt2sq=sigma0(j)**2
2274             sigii1=sigii(i)
2275             sigii2=sigii(j)
2276             ratsig1=sigt2sq/sigt1sq
2277             ratsig2=1.0D0/ratsig1
2278             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2279             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2280             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2281           else
2282             rsum_max=sigma(i,j)
2283           endif
2284 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2285             sigmaii(i,j)=rsum_max
2286             sigmaii(j,i)=rsum_max
2287 !         else
2288 !           sigmaii(i,j)=r0(i,j)
2289 !           sigmaii(j,i)=r0(i,j)
2290 !         endif
2291 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2292           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2293             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2294             augm(i,j)=epsij*r_augm**(2*expon)
2295 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2296             augm(j,i)=augm(i,j)
2297           else
2298             augm(i,j)=0.0D0
2299             augm(j,i)=0.0D0
2300           endif
2301           if (lprint) then
2302             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2303             restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2304             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2305           endif
2306         enddo
2307       enddo
2308
2309       allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2310       allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2311       allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2312       allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2313       allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2314       allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2315       allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2316       allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2317       allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2318       allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2319       allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2320       allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2321       allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2322       allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2323       allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2324       allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2325
2326 !      augm(:,:)=0.0D0
2327 !      chip(:)=0.0D0
2328 !      alp(:)=0.0D0
2329 !      sigma0(:)=0.0D0
2330 !      sigii(:)=0.0D0
2331 !      rr0(:)=0.0D0
2332
2333       read (isidep_nucl,*) ipot_nucl
2334 !      print *,"TU?!",ipot_nucl
2335       if (ipot_nucl.eq.1) then
2336         do i=1,ntyp_molec(2)
2337           do j=i,ntyp_molec(2)
2338             read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2339             elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2340           enddo
2341         enddo
2342       else
2343         do i=1,ntyp_molec(2)
2344           do j=i,ntyp_molec(2)
2345             read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2346                chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2347                elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2348           enddo
2349         enddo
2350       endif
2351 !      rpp(1,1)=2**(1.0/6.0)*5.16158
2352       do i=1,ntyp_molec(2)
2353         do j=i,ntyp_molec(2)
2354           rrij=sigma_nucl(i,j)
2355           r0_nucl(i,j)=rrij
2356           r0_nucl(j,i)=rrij
2357           rrij=rrij**expon
2358           epsij=4*eps_nucl(i,j)
2359           sigeps=dsign(1.0D0,epsij)
2360           epsij=dabs(epsij)
2361           aa_nucl(i,j)=epsij*rrij*rrij
2362           bb_nucl(i,j)=-sigeps*epsij*rrij
2363           ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2364           ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2365           ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2366           ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2367           sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2368          2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2369         enddo
2370         do j=1,i-1
2371           aa_nucl(i,j)=aa_nucl(j,i)
2372           bb_nucl(i,j)=bb_nucl(j,i)
2373           ael3_nucl(i,j)=ael3_nucl(j,i)
2374           ael6_nucl(i,j)=ael6_nucl(j,i)
2375           ael63_nucl(i,j)=ael63_nucl(j,i)
2376           ael32_nucl(i,j)=ael32_nucl(j,i)
2377           elpp3_nucl(i,j)=elpp3_nucl(j,i)
2378           elpp6_nucl(i,j)=elpp6_nucl(j,i)
2379           elpp63_nucl(i,j)=elpp63_nucl(j,i)
2380           elpp32_nucl(i,j)=elpp32_nucl(j,i)
2381           eps_nucl(i,j)=eps_nucl(j,i)
2382           sigma_nucl(i,j)=sigma_nucl(j,i)
2383           sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2384         enddo
2385       enddo
2386
2387       write(iout,*) "tube param"
2388       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2389       ccavtubpep,dcavtubpep,tubetranenepep
2390       sigmapeptube=sigmapeptube**6
2391       sigeps=dsign(1.0D0,epspeptube)
2392       epspeptube=dabs(epspeptube)
2393       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2394       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2395       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2396       do i=1,ntyp
2397        read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2398       ccavtub(i),dcavtub(i),tubetranene(i)
2399        sigmasctube=sigmasctube**6
2400        sigeps=dsign(1.0D0,epssctube)
2401        epssctube=dabs(epssctube)
2402        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2403        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2404       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2405       enddo
2406 !-----------------READING SC BASE POTENTIALS-----------------------------
2407       allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2408       allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2409       allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2410       allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2411       allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2412       allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2413       allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2414       allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2415       allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2416       allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2417       allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2418       allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2419       allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2420       allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2421       allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2422       allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2423
2424
2425       do i=1,ntyp_molec(1)
2426        do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2427         write (*,*) "Im in ", i, " ", j
2428         read(isidep_scbase,*) &
2429         eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2430         chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2431          write(*,*) "eps",eps_scbase(i,j)
2432         read(isidep_scbase,*) &
2433        (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2434        chis_scbase(i,j,1),chis_scbase(i,j,2)
2435         read(isidep_scbase,*) &
2436        dhead_scbasei(i,j), &
2437        dhead_scbasej(i,j), &
2438        rborn_scbasei(i,j),rborn_scbasej(j,i)
2439         read(isidep_scbase,*) &
2440        (wdipdip_scbase(k,i,j),k=1,3), &
2441        (wqdip_scbase(k,i,j),k=1,2)
2442         read(isidep_scbase,*) &
2443        alphapol_scbase(i,j), &
2444        epsintab_scbase(i,j)
2445        END DO
2446       END DO
2447       allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2448       allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2449
2450       do i=1,ntyp_molec(1)
2451        do j=1,ntyp_molec(2)-1
2452           epsij=eps_scbase(i,j)
2453           rrij=sigma_scbase(i,j)
2454 !          r0(i,j)=rrij
2455 !          r0(j,i)=rrij
2456           rrij=rrij**expon
2457 !          epsij=eps(i,j)
2458           sigeps=dsign(1.0D0,epsij)
2459           epsij=dabs(epsij)
2460           aa_scbase(i,j)=epsij*rrij*rrij
2461           bb_scbase(i,j)=-sigeps*epsij*rrij
2462         enddo
2463        enddo
2464 !-----------------READING PEP BASE POTENTIALS-------------------
2465       allocate(eps_pepbase(ntyp_molec(2)))
2466       allocate(sigma_pepbase(ntyp_molec(2)))
2467       allocate(chi_pepbase(ntyp_molec(2),2))
2468       allocate(chipp_pepbase(ntyp_molec(2),2))
2469       allocate(alphasur_pepbase(4,ntyp_molec(2)))
2470       allocate(sigmap1_pepbase(ntyp_molec(2)))
2471       allocate(sigmap2_pepbase(ntyp_molec(2)))
2472       allocate(chis_pepbase(ntyp_molec(2),2))
2473       allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2474
2475
2476        do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2477         write (*,*) "Im in ", i, " ", j
2478         read(isidep_pepbase,*) &
2479         eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2480         chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2481          write(*,*) "eps",eps_pepbase(j)
2482         read(isidep_pepbase,*) &
2483        (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2484        chis_pepbase(j,1),chis_pepbase(j,2)
2485         read(isidep_pepbase,*) &
2486        (wdipdip_pepbase(k,j),k=1,3)
2487        END DO
2488       allocate(aa_pepbase(ntyp_molec(2)))
2489       allocate(bb_pepbase(ntyp_molec(2)))
2490
2491        do j=1,ntyp_molec(2)-1
2492           epsij=eps_pepbase(j)
2493           rrij=sigma_pepbase(j)
2494 !          r0(i,j)=rrij
2495 !          r0(j,i)=rrij
2496           rrij=rrij**expon
2497 !          epsij=eps(i,j)
2498           sigeps=dsign(1.0D0,epsij)
2499           epsij=dabs(epsij)
2500           aa_pepbase(j)=epsij*rrij*rrij
2501           bb_pepbase(j)=-sigeps*epsij*rrij
2502         enddo
2503 !--------------READING SC PHOSPHATE-------------------------------------
2504       allocate(eps_scpho(ntyp_molec(1)))
2505       allocate(sigma_scpho(ntyp_molec(1)))
2506       allocate(chi_scpho(ntyp_molec(1),2))
2507       allocate(chipp_scpho(ntyp_molec(1),2))
2508       allocate(alphasur_scpho(4,ntyp_molec(1)))
2509       allocate(sigmap1_scpho(ntyp_molec(1)))
2510       allocate(sigmap2_scpho(ntyp_molec(1)))
2511       allocate(chis_scpho(ntyp_molec(1),2))
2512       allocate(wqq_scpho(ntyp_molec(1)))
2513       allocate(wqdip_scpho(2,ntyp_molec(1)))
2514       allocate(alphapol_scpho(ntyp_molec(1)))
2515       allocate(epsintab_scpho(ntyp_molec(1)))
2516       allocate(dhead_scphoi(ntyp_molec(1)))
2517       allocate(rborn_scphoi(ntyp_molec(1)))
2518       allocate(rborn_scphoj(ntyp_molec(1)))
2519
2520 !      j=1
2521        do j=1,ntyp_molec(1) ! without U then we will take T for U
2522         write (*,*) "Im in scpho ", i, " ", j
2523         read(isidep_scpho,*) &
2524         eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2525         chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2526          write(*,*) "eps",eps_scpho(j)
2527         read(isidep_scpho,*) &
2528        (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2529        chis_scpho(j,1),chis_scpho(j,2)
2530         read(isidep_scpho,*) &
2531        (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2532         read(isidep_scpho,*) &
2533          epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j)
2534
2535        END DO
2536       allocate(aa_scpho(ntyp_molec(1)))
2537       allocate(bb_scpho(ntyp_molec(1)))
2538
2539        do j=1,ntyp_molec(1)
2540           epsij=eps_scpho(j)
2541           rrij=sigma_scpho(j)
2542 !          r0(i,j)=rrij
2543 !          r0(j,i)=rrij
2544           rrij=rrij**expon
2545 !          epsij=eps(i,j)
2546           sigeps=dsign(1.0D0,epsij)
2547           epsij=dabs(epsij)
2548           aa_scpho(j)=epsij*rrij*rrij
2549           bb_scpho(j)=-sigeps*epsij*rrij
2550         enddo
2551
2552
2553         read(isidep_peppho,*) &
2554         eps_peppho,sigma_peppho
2555         read(isidep_peppho,*) &
2556        (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2557         read(isidep_peppho,*) &
2558        (wqdip_peppho(k),k=1,2)
2559
2560           epsij=eps_peppho
2561           rrij=sigma_peppho
2562 !          r0(i,j)=rrij
2563 !          r0(j,i)=rrij
2564           rrij=rrij**expon
2565 !          epsij=eps(i,j)
2566           sigeps=dsign(1.0D0,epsij)
2567           epsij=dabs(epsij)
2568           aa_peppho=epsij*rrij*rrij
2569           bb_peppho=-sigeps*epsij*rrij
2570
2571
2572       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2573       bad(:,:)=0.0D0
2574
2575 #ifdef OLDSCP
2576 !
2577 ! Define the SC-p interaction constants (hard-coded; old style)
2578 !
2579       do i=1,ntyp
2580 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2581 ! helix formation)
2582 !       aad(i,1)=0.3D0*4.0D0**12
2583 ! Following line for constants currently implemented
2584 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2585         aad(i,1)=1.5D0*4.0D0**12
2586 !       aad(i,1)=0.17D0*5.6D0**12
2587         aad(i,2)=aad(i,1)
2588 ! "Soft" SC-p repulsion
2589         bad(i,1)=0.0D0
2590 ! Following line for constants currently implemented
2591 !       aad(i,1)=0.3D0*4.0D0**6
2592 ! "Hard" SC-p repulsion
2593         bad(i,1)=3.0D0*4.0D0**6
2594 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2595         bad(i,2)=bad(i,1)
2596 !       aad(i,1)=0.0D0
2597 !       aad(i,2)=0.0D0
2598 !       bad(i,1)=1228.8D0
2599 !       bad(i,2)=1228.8D0
2600       enddo
2601 #else
2602 !
2603 ! 8/9/01 Read the SC-p interaction constants from file
2604 !
2605       do i=1,ntyp
2606         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2607       enddo
2608       do i=1,ntyp
2609         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2610         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2611         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2612         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2613       enddo
2614 !      lprint=.true.
2615       if (lprint) then
2616         write (iout,*) "Parameters of SC-p interactions:"
2617         do i=1,ntyp
2618           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2619            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2620         enddo
2621       endif
2622 !      lprint=.false.
2623 #endif
2624       allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2625
2626       do i=1,ntyp_molec(2)
2627         read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2628       enddo
2629       do i=1,ntyp_molec(2)
2630         aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2631         bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2632       enddo
2633       r0pp=1.12246204830937298142*5.16158
2634       epspp=4.95713/4
2635       AEES=108.661
2636       BEES=0.433246
2637
2638 !
2639 ! Define the constants of the disulfide bridge
2640 !
2641       ebr=-5.50D0
2642 !
2643 ! Old arbitrary potential - commented out.
2644 !
2645 !      dbr= 4.20D0
2646 !      fbr= 3.30D0
2647 !
2648 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2649 ! energy surface of diethyl disulfide.
2650 ! A. Liwo and U. Kozlowska, 11/24/03
2651 !
2652       D0CM = 3.78d0
2653       AKCM = 15.1d0
2654       AKTH = 11.0d0
2655       AKCT = 12.0d0
2656       V1SS =-1.08d0
2657       V2SS = 7.61d0
2658       V3SS = 13.7d0
2659 !      akcm=0.0d0
2660 !      akth=0.0d0
2661 !      akct=0.0d0
2662 !      v1ss=0.0d0
2663 !      v2ss=0.0d0
2664 !      v3ss=0.0d0
2665
2666       if(me.eq.king) then
2667       write (iout,'(/a)') "Disulfide bridge parameters:"
2668       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2669       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2670       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2671       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2672         ' v3ss:',v3ss
2673       endif
2674       return
2675   111 write (iout,*) "Error reading bending energy parameters."
2676       goto 999
2677   112 write (iout,*) "Error reading rotamer energy parameters."
2678       goto 999
2679   113 write (iout,*) "Error reading torsional energy parameters."
2680       goto 999
2681   114 write (iout,*) "Error reading double torsional energy parameters."
2682       goto 999
2683   115 write (iout,*) &
2684         "Error reading cumulant (multibody energy) parameters."
2685       goto 999
2686   116 write (iout,*) "Error reading electrostatic energy parameters."
2687       goto 999
2688   117 write (iout,*) "Error reading side chain interaction parameters."
2689       goto 999
2690   118 write (iout,*) "Error reading SCp interaction parameters."
2691       goto 999
2692   119 write (iout,*) "Error reading SCCOR parameters"
2693   999 continue
2694 #ifdef MPI
2695       call MPI_Finalize(Ierror)
2696 #endif
2697       stop
2698       return
2699       end subroutine parmread
2700 #endif
2701 !-----------------------------------------------------------------------------
2702 ! printmat.f
2703 !-----------------------------------------------------------------------------
2704       subroutine printmat(ldim,m,n,iout,key,a)
2705
2706       integer :: n,ldim
2707       character(len=3),dimension(n) :: key
2708       real(kind=8),dimension(ldim,n) :: a
2709 !el local variables
2710       integer :: i,j,k,m,iout,nlim
2711
2712       do 1 i=1,n,8
2713       nlim=min0(i+7,n)
2714       write (iout,1000) (key(k),k=i,nlim)
2715       write (iout,1020)
2716  1000 format (/5x,8(6x,a3))
2717  1020 format (/80(1h-)/)
2718       do 2 j=1,n
2719       write (iout,1010) key(j),(a(j,k),k=i,nlim)
2720     2 continue
2721     1 continue
2722  1010 format (a3,2x,8(f9.4))
2723       return
2724       end subroutine printmat
2725 !-----------------------------------------------------------------------------
2726 ! readpdb.F
2727 !-----------------------------------------------------------------------------
2728       subroutine readpdb
2729 ! Read the PDB file and convert the peptide geometry into virtual-chain
2730 ! geometry.
2731       use geometry_data
2732       use energy_data, only: itype,istype
2733       use control_data
2734       use compare_data
2735       use MPI_data
2736       use control, only: rescode,sugarcode
2737 !      implicit real*8 (a-h,o-z)
2738 !      include 'DIMENSIONS'
2739 !      include 'COMMON.LOCAL'
2740 !      include 'COMMON.VAR'
2741 !      include 'COMMON.CHAIN'
2742 !      include 'COMMON.INTERACT'
2743 !      include 'COMMON.IOUNITS'
2744 !      include 'COMMON.GEO'
2745 !      include 'COMMON.NAMES'
2746 !      include 'COMMON.CONTROL'
2747 !      include 'COMMON.DISTFIT'
2748 !      include 'COMMON.SETUP'
2749       integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2750 !        ishift_pdb
2751       logical :: lprn=.true.,fail
2752       real(kind=8),dimension(3) :: e1,e2,e3
2753       real(kind=8) :: dcj,efree_temp
2754       character(len=3) :: seq,res
2755       character(len=5) :: atom
2756       character(len=80) :: card
2757       real(kind=8),dimension(3,20) :: sccor
2758       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin       !rescode,
2759       integer :: isugar,molecprev,firstion
2760       character*1 :: sugar
2761       real(kind=8) :: cou
2762       real(kind=8),dimension(3) :: ccc
2763 !el local varables
2764       integer,dimension(2,maxres/3) :: hfrag_alloc
2765       integer,dimension(4,maxres/3) :: bfrag_alloc
2766       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2767       real(kind=8),dimension(:,:), allocatable  :: c_temporary
2768       integer,dimension(:,:) , allocatable  :: itype_temporary
2769       integer,dimension(:),allocatable :: istype_temp
2770       efree_temp=0.0d0
2771       ibeg=1
2772       ishift1=0
2773       ishift=0
2774       molecule=1
2775       counter=0
2776 !      write (2,*) "UNRES_PDB",unres_pdb
2777       ires=0
2778       ires_old=0
2779       nres=0
2780       iii=0
2781       lsecondary=.false.
2782       nhfrag=0
2783       nbfrag=0
2784 !-----------------------------
2785       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2786       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2787
2788       do i=1,100000
2789         read (ipdbin,'(a80)',end=10) card
2790 !       write (iout,'(a)') card
2791         if (card(:5).eq.'HELIX') then
2792           nhfrag=nhfrag+1
2793           lsecondary=.true.
2794           read(card(22:25),*) hfrag(1,nhfrag)
2795           read(card(34:37),*) hfrag(2,nhfrag)
2796         endif
2797         if (card(:5).eq.'SHEET') then
2798           nbfrag=nbfrag+1
2799           lsecondary=.true.
2800           read(card(24:26),*) bfrag(1,nbfrag)
2801           read(card(35:37),*) bfrag(2,nbfrag)
2802 !rc----------------------------------------
2803 !rc  to be corrected !!!
2804           bfrag(3,nbfrag)=bfrag(1,nbfrag)
2805           bfrag(4,nbfrag)=bfrag(2,nbfrag)
2806 !rc----------------------------------------
2807         endif
2808         if (card(:3).eq.'END') then
2809           goto 10
2810         else if (card(:3).eq.'TER') then
2811 ! End current chain
2812           ires_old=ires+2
2813           ishift=ishift+1
2814           ishift1=ishift1+1
2815           itype(ires_old,molecule)=ntyp1_molec(molecule)
2816           itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2817           nres_molec(molecule)=nres_molec(molecule)+2
2818           ibeg=2
2819 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2820           if (unres_pdb) then
2821             do j=1,3
2822               dc(j,ires)=sccor(j,iii)
2823             enddo
2824           else
2825             call sccenter(ires,iii,sccor)
2826 !          iii=0
2827           endif
2828           iii=0
2829         endif
2830 ! Read free energy
2831         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2832 ! Fish out the ATOM cards.
2833 !        write(iout,*) 'card',card(1:20)
2834         if (index(card(1:4),'ATOM').gt.0) then
2835           read (card(12:16),*) atom
2836 !          write (iout,*) "! ",atom," !",ires
2837 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2838           read (card(23:26),*) ires
2839           read (card(18:20),'(a3)') res
2840 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2841 !     &      " ires_old",ires_old
2842 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2843 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2844           if (ires-ishift+ishift1.ne.ires_old) then
2845 ! Calculate the CM of the preceding residue.
2846 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2847             if (ibeg.eq.0) then
2848 !              write (iout,*) "Calculating sidechain center iii",iii
2849               if (unres_pdb) then
2850                 do j=1,3
2851                   dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2852                 enddo
2853               else
2854                 call sccenter(ires_old,iii,sccor)
2855               endif
2856               iii=0
2857             endif
2858 ! Start new residue.
2859             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2860               ires=ires_old
2861               cycle
2862             else if (ibeg.eq.1) then
2863               write (iout,*) "BEG ires",ires
2864               ishift=ires-1
2865               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2866                 ishift=ishift-1
2867                 itype(1,1)=ntyp1
2868                 nres_molec(molecule)=nres_molec(molecule)+1
2869               endif
2870               ires=ires-ishift+ishift1
2871               ires_old=ires
2872 !              write (iout,*) "ishift",ishift," ires",ires,&
2873 !               " ires_old",ires_old
2874               ibeg=0
2875             else if (ibeg.eq.2) then
2876 ! Start a new chain
2877               ishift=-ires_old+ires-1 !!!!!
2878               ishift1=ishift1-1    !!!!!
2879 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2880               ires=ires-ishift+ishift1
2881 !              print *,ires,ishift,ishift1
2882               ires_old=ires
2883               ibeg=0
2884             else
2885               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2886               ires=ires-ishift+ishift1
2887               ires_old=ires
2888             endif
2889 !            print *,'atom',ires,atom
2890             if (res.eq.'ACE' .or. res.eq.'NHE') then
2891               itype(ires,1)=10
2892             else
2893              if (atom.eq.'CA  '.or.atom.eq.'N   ') then
2894              molecule=1
2895               itype(ires,molecule)=rescode(ires,res,0,molecule)
2896               firstion=0
2897 !              nres_molec(molecule)=nres_molec(molecule)+1
2898             else
2899              molecule=2
2900               itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2901 !              nres_molec(molecule)=nres_molec(molecule)+1
2902              read (card(19:19),'(a1)') sugar
2903              isugar=sugarcode(sugar,ires)
2904 !            if (ibeg.eq.1) then
2905 !              istype(1)=isugar
2906 !            else
2907               istype(ires)=isugar
2908 !              print *,"ires=",ires,istype(ires)
2909 !            endif
2910
2911             endif
2912             endif
2913           else
2914             ires=ires-ishift+ishift1
2915           endif
2916 !          write (iout,*) "ires_old",ires_old," ires",ires
2917           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2918 !            ishift1=ishift1+1
2919           endif
2920 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2921           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2922              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2923             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2924 !              print *,ires,ishift,ishift1
2925 !            write (iout,*) "backbone ",atom
2926 #ifdef DEBUG
2927             write (iout,'(2i3,2x,a,3f8.3)') &
2928             ires,itype(ires,1),res,(c(j,ires),j=1,3)
2929 #endif
2930             iii=iii+1
2931               nres_molec(molecule)=nres_molec(molecule)+1
2932             do j=1,3
2933               sccor(j,iii)=c(j,ires)
2934             enddo
2935           else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2936                atom.eq."C2'" .or. atom.eq."C3'" &
2937                .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2938             read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2939 !c            write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2940 !              print *,ires,ishift,ishift1
2941             counter=counter+1
2942 !            iii=iii+1
2943 !            do j=1,3
2944 !              sccor(j,iii)=c(j,ires)
2945 !            enddo
2946             do j=1,3
2947               c(j,ires)=c(j,ires)+ccc(j)/5.0
2948             enddo
2949              print *,counter,molecule
2950              if (counter.eq.5) then
2951 !            iii=iii+1
2952               nres_molec(molecule)=nres_molec(molecule)+1
2953               firstion=0
2954 !            do j=1,3
2955 !              sccor(j,iii)=c(j,ires)
2956 !            enddo
2957              counter=0
2958            endif
2959 !            print *, "ATOM",atom(1:3)
2960           else if (atom.eq."C5'") then
2961              read (card(19:19),'(a1)') sugar
2962              isugar=sugarcode(sugar,ires)
2963             if (ibeg.eq.1) then
2964               istype(1)=isugar
2965             else
2966               istype(ires)=isugar
2967 !              print *,ires,istype(ires)
2968             endif
2969             if (unres_pdb) then
2970               read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2971             else
2972               iii=iii+1
2973               read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2974             endif
2975 !            write (*,*) card(23:27),ires,itype(ires,1)
2976           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2977                    atom.ne.'N' .and. atom.ne.'C' .and. &
2978                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2979                    atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2980                    .and. atom.ne.'P  '.and. &
2981                   atom(1:1).ne.'H' .and. &
2982                   atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2983                   ) then
2984 !            write (iout,*) "sidechain ",atom
2985 !            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2986                  if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2987 !                        write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2988             iii=iii+1
2989             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2990               endif
2991           endif
2992         else if ((ions).and.(card(1:6).eq.'HETATM')) then
2993        if (firstion.eq.0) then
2994        firstion=1
2995        if (unres_pdb) then
2996          do j=1,3
2997            dc(j,ires)=sccor(j,iii)
2998          enddo
2999        else
3000           call sccenter(ires,iii,sccor)
3001        endif
3002        endif
3003           read (card(12:16),*) atom
3004           print *,"HETATOM", atom
3005           read (card(18:20),'(a3)') res
3006           if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3007           (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
3008           .or.(atom(1:2).eq.'K ')) &
3009           then
3010            ires=ires+1
3011            if (molecule.ne.5) molecprev=molecule
3012            molecule=5
3013            nres_molec(molecule)=nres_molec(molecule)+1
3014            print *,"HERE",nres_molec(molecule)
3015            itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
3016            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3017           endif
3018         endif !atom
3019       enddo
3020    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3021       if (ires.eq.0) return
3022 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3023 ! system
3024       nres=ires
3025       if ((ires_old.ne.ires).and.(molecule.ne.5)) &
3026          nres_molec(molecule)=nres_molec(molecule)-2
3027 !      print *,'I have', nres_molec(:)
3028
3029       do k=1,4 ! ions are without dummy
3030        if (nres_molec(k).eq.0) cycle
3031       do i=2,nres-1
3032 !        write (iout,*) i,itype(i,1)
3033 !        if (itype(i,1).eq.ntyp1) then
3034 !          write (iout,*) "dummy",i,itype(i,1)
3035 !          do j=1,3
3036 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3037 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
3038 !            dc(j,i)=c(j,i)
3039 !          enddo
3040 !        endif
3041         if (itype(i,k).eq.ntyp1_molec(k)) then
3042          if (itype(i+1,k).eq.ntyp1_molec(k)) then
3043           if (itype(i+2,k).eq.0) then
3044 !           print *,"masz sieczke"
3045            do j=1,5
3046             if (itype(i+2,j).ne.0) then
3047             itype(i+1,k)=0
3048             itype(i+1,j)=ntyp1_molec(j)
3049             nres_molec(k)=nres_molec(k)-1
3050             nres_molec(j)=nres_molec(j)+1
3051             go to 3331
3052             endif
3053            enddo
3054  3331      continue
3055           endif
3056 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3057 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3058 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3059 !           if (unres_pdb) then
3060 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3061 !            print *,i,'tu dochodze'
3062 !            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3063 !            if (fail) then
3064 !              e2(1)=0.0d0
3065 !              e2(2)=1.0d0
3066 !              e2(3)=0.0d0
3067 !            endif !fail
3068 !            print *,i,'a tu?'
3069 !            do j=1,3
3070 !             c(j,i)=c(j,i-1)-1.9d0*e2(j)
3071 !            enddo
3072 !           else   !unres_pdb
3073            do j=1,3
3074              dcj=(c(j,i-2)-c(j,i-3))/2.0
3075             if (dcj.eq.0) dcj=1.23591524223
3076              c(j,i)=c(j,i-1)+dcj
3077              c(j,nres+i)=c(j,i)
3078            enddo
3079 !          endif   !unres_pdb
3080          else     !itype(i+1,1).eq.ntyp1
3081 !          if (unres_pdb) then
3082 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3083 !            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3084 !            if (fail) then
3085 !              e2(1)=0.0d0
3086 !              e2(2)=1.0d0
3087 !              e2(3)=0.0d0
3088 !            endif
3089             do j=1,3
3090               c(j,i)=c(j,i+1)-1.9d0*e2(j)
3091             enddo
3092 !          else !unres_pdb
3093            do j=1,3
3094             dcj=(c(j,i+3)-c(j,i+2))/2.0
3095             if (dcj.eq.0) dcj=1.23591524223
3096             c(j,i)=c(j,i+1)-dcj
3097             c(j,nres+i)=c(j,i)
3098            enddo
3099 !          endif !unres_pdb
3100          endif !itype(i+1,1).eq.ntyp1
3101         endif  !itype.eq.ntyp1
3102
3103       enddo
3104       enddo
3105 ! Calculate the CM of the last side chain.
3106       if (iii.gt.0)  then
3107       if (unres_pdb) then
3108         do j=1,3
3109           dc(j,ires)=sccor(j,iii)
3110         enddo
3111       else
3112         call sccenter(ires,iii,sccor)
3113       endif
3114       endif
3115 !      nres=ires
3116       nsup=nres
3117       nstart_sup=1
3118 !      print *,"molecule",molecule
3119       if ((itype(nres,1).ne.10)) then
3120         nres=nres+1
3121           if (molecule.eq.5) molecule=molecprev
3122         itype(nres,molecule)=ntyp1_molec(molecule)
3123         nres_molec(molecule)=nres_molec(molecule)+1
3124 !        if (unres_pdb) then
3125 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3126 !          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3127 !          if (fail) then
3128 !            e2(1)=0.0d0
3129 !            e2(2)=1.0d0
3130 !            e2(3)=0.0d0
3131 !          endif
3132 !          do j=1,3
3133 !            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3134 !          enddo
3135 !        else
3136         do j=1,3
3137           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3138           c(j,nres)=c(j,nres-1)+dcj
3139           c(j,2*nres)=c(j,nres)
3140         enddo
3141 !        endif
3142       endif
3143 !     print *,'I have',nres, nres_molec(:)
3144
3145 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3146 #ifdef WHAM_RUN
3147       if (nres.ne.nres0) then
3148         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3149                        " NRES0=",nres0
3150         stop "Error nres value in WHAM input"
3151       endif
3152 #endif
3153 !---------------------------------
3154 !el reallocate tables
3155 !      do i=1,maxres/3
3156 !       do j=1,2
3157 !         hfrag_alloc(j,i)=hfrag(j,i)
3158 !        enddo
3159 !       do j=1,4
3160 !         bfrag_alloc(j,i)=bfrag(j,i)
3161 !        enddo
3162 !      enddo
3163
3164 !      deallocate(hfrag)
3165 !      deallocate(bfrag)
3166 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
3167 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3168 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3169 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
3170
3171 !      do i=1,nhfrag
3172 !       do j=1,2
3173 !         hfrag(j,i)=hfrag_alloc(j,i)
3174 !        enddo
3175 !      enddo
3176 !      do i=1,nbfrag
3177 !       do j=1,4
3178 !         bfrag(j,i)=bfrag_alloc(j,i)
3179 !        enddo
3180 !      enddo
3181 !el end reallocate tables
3182 !---------------------------------
3183       do i=2,nres-1
3184         do j=1,3
3185           c(j,i+nres)=dc(j,i)
3186         enddo
3187       enddo
3188       do j=1,3
3189         c(j,nres+1)=c(j,1)
3190         c(j,2*nres)=c(j,nres)
3191       enddo
3192
3193       if (itype(1,1).eq.ntyp1) then
3194         nsup=nsup-1
3195         nstart_sup=2
3196         if (unres_pdb) then
3197 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3198           call refsys(2,3,4,e1,e2,e3,fail)
3199           if (fail) then
3200             e2(1)=0.0d0
3201             e2(2)=1.0d0
3202             e2(3)=0.0d0
3203           endif
3204           do j=1,3
3205             c(j,1)=c(j,2)-1.9d0*e2(j)
3206           enddo
3207         else
3208         do j=1,3
3209           dcj=(c(j,4)-c(j,3))/2.0
3210           c(j,1)=c(j,2)-dcj
3211           c(j,nres+1)=c(j,1)
3212         enddo
3213         endif
3214       endif
3215 ! First lets assign correct dummy to correct type of chain
3216 ! 1) First residue
3217       if (itype(1,1).eq.ntyp1) then
3218         if (itype(2,1).eq.0) then
3219          do j=2,5
3220            if (itype(2,j).ne.0) then
3221            itype(1,1)=0
3222            itype(1,j)=ntyp1_molec(j)
3223            nres_molec(1)=nres_molec(1)-1
3224            nres_molec(j)=nres_molec(j)+1
3225            go to 3231
3226            endif
3227          enddo
3228 3231    continue
3229         endif
3230        endif
3231        print *,'I have',nres, nres_molec(:)
3232
3233 ! Copy the coordinates to reference coordinates
3234 !      do i=1,2*nres
3235 !        do j=1,3
3236 !          cref(j,i)=c(j,i)
3237 !        enddo
3238 !      enddo
3239 ! Calculate internal coordinates.
3240       if (lprn) then
3241       write (iout,'(/a)') &
3242         "Cartesian coordinates of the reference structure"
3243       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3244        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3245       do ires=1,nres
3246         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3247           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3248           (c(j,ires+nres),j=1,3)
3249       enddo
3250       endif
3251 ! znamy już nres wiec mozna alokowac tablice
3252 ! Calculate internal coordinates.
3253       if(me.eq.king.or..not.out1file)then
3254        write (iout,'(a)') &
3255          "Backbone and SC coordinates as read from the PDB"
3256        do ires=1,nres
3257         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3258           ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3259           (c(j,nres+ires),j=1,3)
3260        enddo
3261       endif
3262 ! NOW LETS ROCK! SORTING
3263       allocate(c_temporary(3,2*nres))
3264       allocate(itype_temporary(nres,5))
3265       allocate(molnum(nres+1))
3266       allocate(istype_temp(nres))
3267        itype_temporary(:,:)=0
3268       seqalingbegin=1
3269       do k=1,5
3270         do i=1,nres
3271          if (itype(i,k).ne.0) then
3272           do j=1,3
3273           c_temporary(j,seqalingbegin)=c(j,i)
3274           c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3275
3276           enddo
3277           itype_temporary(seqalingbegin,k)=itype(i,k)
3278           print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3279           istype_temp(seqalingbegin)=istype(i)
3280           molnum(seqalingbegin)=k
3281           seqalingbegin=seqalingbegin+1
3282          endif
3283         enddo
3284        enddo
3285        do i=1,2*nres
3286         do j=1,3
3287         c(j,i)=c_temporary(j,i)
3288         enddo
3289        enddo
3290        do k=1,5
3291         do i=1,nres
3292          itype(i,k)=itype_temporary(i,k)
3293          istype(i)=istype_temp(i)
3294         enddo
3295        enddo
3296       if (itype(1,1).eq.ntyp1) then
3297         nsup=nsup-1
3298         nstart_sup=2
3299         if (unres_pdb) then
3300 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3301           call refsys(2,3,4,e1,e2,e3,fail)
3302           if (fail) then
3303             e2(1)=0.0d0
3304             e2(2)=1.0d0
3305             e2(3)=0.0d0
3306           endif
3307           do j=1,3
3308             c(j,1)=c(j,2)-1.9d0*e2(j)
3309           enddo
3310         else
3311         do j=1,3
3312           dcj=(c(j,4)-c(j,3))/2.0
3313           c(j,1)=c(j,2)-dcj
3314           c(j,nres+1)=c(j,1)
3315         enddo
3316         endif
3317       endif
3318
3319       if (lprn) then
3320       write (iout,'(/a)') &
3321         "Cartesian coordinates of the reference structure after sorting"
3322       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3323        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3324       do ires=1,nres
3325         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3326           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3327           (c(j,ires+nres),j=1,3)
3328       enddo
3329       endif
3330
3331 !       print *,seqalingbegin,nres
3332       if(.not.allocated(vbld)) then
3333        allocate(vbld(2*nres))
3334        do i=1,2*nres
3335          vbld(i)=0.d0
3336        enddo
3337       endif
3338       if(.not.allocated(vbld_inv)) then
3339        allocate(vbld_inv(2*nres))
3340        do i=1,2*nres
3341          vbld_inv(i)=0.d0
3342        enddo
3343       endif
3344 !!!el
3345       if(.not.allocated(theta)) then
3346         allocate(theta(nres+2))
3347         theta(:)=0.0d0
3348       endif
3349
3350       if(.not.allocated(phi)) allocate(phi(nres+2))
3351       if(.not.allocated(alph)) allocate(alph(nres+2))
3352       if(.not.allocated(omeg)) allocate(omeg(nres+2))
3353       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3354       if(.not.allocated(phiref)) allocate(phiref(nres+2))
3355       if(.not.allocated(costtab)) allocate(costtab(nres))
3356       if(.not.allocated(sinttab)) allocate(sinttab(nres))
3357       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3358       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3359       if(.not.allocated(xxref)) allocate(xxref(nres))
3360       if(.not.allocated(yyref)) allocate(yyref(nres))
3361       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3362       if(.not.allocated(dc_norm)) then
3363 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3364         allocate(dc_norm(3,0:2*nres+2))
3365         dc_norm(:,:)=0.d0
3366       endif
3367
3368       call int_from_cart(.true.,.false.)
3369       call sc_loc_geom(.false.)
3370       do i=1,nres
3371         thetaref(i)=theta(i)
3372         phiref(i)=phi(i)
3373       enddo
3374 !      do i=1,2*nres
3375 !        vbld_inv(i)=0.d0
3376 !        vbld(i)=0.d0
3377 !      enddo
3378
3379       do i=1,nres-1
3380         do j=1,3
3381           dc(j,i)=c(j,i+1)-c(j,i)
3382           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3383         enddo
3384       enddo
3385       do i=2,nres-1
3386         do j=1,3
3387           dc(j,i+nres)=c(j,i+nres)-c(j,i)
3388           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3389         enddo
3390 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3391 !        vbld_inv(i+nres)
3392       enddo
3393 !      call chainbuild
3394 ! Copy the coordinates to reference coordinates
3395 ! Splits to single chain if occurs
3396
3397 !      do i=1,2*nres
3398 !        do j=1,3
3399 !          cref(j,i,cou)=c(j,i)
3400 !        enddo
3401 !      enddo
3402 !
3403       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3404       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3405       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3406 !-----------------------------
3407       kkk=1
3408       lll=0
3409       cou=1
3410         write (iout,*) "symetr", symetr
3411       do i=1,nres
3412       lll=lll+1
3413 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3414       if (i.gt.1) then
3415       if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3416       chain_length=lll-1
3417       kkk=kkk+1
3418 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3419       lll=1
3420       endif
3421       endif
3422         do j=1,3
3423           cref(j,i,cou)=c(j,i)
3424           cref(j,i+nres,cou)=c(j,i+nres)
3425           if (i.le.nres) then
3426           chain_rep(j,lll,kkk)=c(j,i)
3427           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3428           endif
3429          enddo
3430       enddo
3431       write (iout,*) chain_length
3432       if (chain_length.eq.0) chain_length=nres
3433       do j=1,3
3434       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3435       chain_rep(j,chain_length+nres,symetr) &
3436       =chain_rep(j,chain_length+nres,1)
3437       enddo
3438 ! diagnostic
3439 !       write (iout,*) "spraw lancuchy",chain_length,symetr
3440 !       do i=1,4
3441 !         do kkk=1,chain_length
3442 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3443 !         enddo
3444 !        enddo
3445 ! enddiagnostic
3446 ! makes copy of chains
3447         write (iout,*) "symetr", symetr
3448       do j=1,3
3449       dc(j,0)=c(j,1)
3450       enddo
3451
3452       if (symetr.gt.1) then
3453        call permut(symetr)
3454        nperm=1
3455        do i=1,symetr
3456        nperm=nperm*i
3457        enddo
3458        do i=1,nperm
3459        write(iout,*) (tabperm(i,kkk),kkk=1,4)
3460        enddo
3461        do i=1,nperm
3462         cou=0
3463         do kkk=1,symetr
3464          icha=tabperm(i,kkk)
3465 !         write (iout,*) i,icha
3466          do lll=1,chain_length
3467           cou=cou+1
3468            if (cou.le.nres) then
3469            do j=1,3
3470             kupa=mod(lll,chain_length)
3471             iprzes=(kkk-1)*chain_length+lll
3472             if (kupa.eq.0) kupa=chain_length
3473 !            write (iout,*) "kupa", kupa
3474             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3475             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3476           enddo
3477           endif
3478          enddo
3479         enddo
3480        enddo
3481        endif
3482 !-koniec robienia kopii
3483 ! diag
3484       do kkk=1,nperm
3485       write (iout,*) "nowa struktura", nperm
3486       do i=1,nres
3487         write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3488       cref(2,i,kkk),&
3489       cref(3,i,kkk),cref(1,nres+i,kkk),&
3490       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3491       enddo
3492   100 format (//'              alpha-carbon coordinates       ',&
3493                 '     centroid coordinates'/ &
3494                 '       ', 6X,'X',11X,'Y',11X,'Z', &
3495                                 10X,'X',11X,'Y',11X,'Z')
3496   110 format (a,'(',i3,')',6f12.5)
3497
3498       enddo
3499 !c enddiag
3500       do j=1,nbfrag
3501         do i=1,4
3502          bfrag(i,j)=bfrag(i,j)-ishift
3503         enddo
3504       enddo
3505
3506       do j=1,nhfrag
3507         do i=1,2
3508          hfrag(i,j)=hfrag(i,j)-ishift
3509         enddo
3510       enddo
3511       ishift_pdb=ishift
3512
3513       return
3514       end subroutine readpdb
3515 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3516 !-----------------------------------------------------------------------------
3517 ! readrtns_CSA.F
3518 !-----------------------------------------------------------------------------
3519       subroutine read_control
3520 !
3521 ! Read contorl data
3522 !
3523 !      use geometry_data
3524       use comm_machsw
3525       use energy_data
3526       use control_data
3527       use compare_data
3528       use MCM_data
3529       use map_data
3530       use csa_data
3531       use MD_data
3532       use MPI_data
3533       use random, only: random_init
3534 !      implicit real*8 (a-h,o-z)
3535 !      include 'DIMENSIONS'
3536 #ifdef MP
3537       use prng, only:prng_restart
3538       include 'mpif.h'
3539       logical :: OKRandom!, prng_restart
3540       real(kind=8) :: r1
3541 #endif
3542 !      include 'COMMON.IOUNITS'
3543 !      include 'COMMON.TIME1'
3544 !      include 'COMMON.THREAD'
3545 !      include 'COMMON.SBRIDGE'
3546 !      include 'COMMON.CONTROL'
3547 !      include 'COMMON.MCM'
3548 !      include 'COMMON.MAP'
3549 !      include 'COMMON.HEADER'
3550 !      include 'COMMON.CSA'
3551 !      include 'COMMON.CHAIN'
3552 !      include 'COMMON.MUCA'
3553 !      include 'COMMON.MD'
3554 !      include 'COMMON.FFIELD'
3555 !      include 'COMMON.INTERACT'
3556 !      include 'COMMON.SETUP'
3557 !el      integer :: KDIAG,ICORFL,IXDR
3558 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3559       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3560         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
3561 !      character(len=80) :: ucase
3562       character(len=640) :: controlcard
3563
3564       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3565       integer i
3566
3567       nglob_csa=0
3568       eglob_csa=1d99
3569       nmin_csa=0
3570       read (INP,'(a)') titel
3571       call card_concat(controlcard,.true.)
3572 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3573 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3574       call reada(controlcard,'SEED',seed,0.0D0)
3575       call random_init(seed)
3576 ! Set up the time limit (caution! The time must be input in minutes!)
3577       read_cart=index(controlcard,'READ_CART').gt.0
3578       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3579       call readi(controlcard,'SYM',symetr,1)
3580       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3581       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3582       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3583       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3584       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3585       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3586       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3587       call reada(controlcard,'DRMS',drms,0.1D0)
3588       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3589        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3590        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3591        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3592        write (iout,'(a,f10.1)')'DRMS    = ',drms
3593        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3594        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3595       endif
3596       call readi(controlcard,'NZ_START',nz_start,0)
3597       call readi(controlcard,'NZ_END',nz_end,0)
3598       call readi(controlcard,'IZ_SC',iz_sc,0)
3599       timlim=60.0D0*timlim
3600       safety = 60.0d0*safety
3601       timem=timlim
3602       modecalc=0
3603       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3604 !C SHIELD keyword sets if the shielding effect of side-chains is used
3605 !C 0 denots no shielding is used all peptide are equally despite the
3606 !C solvent accesible area
3607 !C 1 the newly introduced function
3608 !C 2 reseved for further possible developement
3609       call readi(controlcard,'SHIELD',shield_mode,0)
3610 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3611         write(iout,*) "shield_mode",shield_mode
3612 !C  Varibles set size of box
3613       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3614       protein=index(controlcard,"PROTEIN").gt.0
3615       ions=index(controlcard,"IONS").gt.0
3616       nucleic=index(controlcard,"NUCLEIC").gt.0
3617       write (iout,*) "with_theta_constr ",with_theta_constr
3618       AFMlog=(index(controlcard,'AFM'))
3619       selfguide=(index(controlcard,'SELFGUIDE'))
3620       print *,'AFMlog',AFMlog,selfguide,"KUPA"
3621       call readi(controlcard,'GENCONSTR',genconstr,0)
3622       call reada(controlcard,'BOXX',boxxsize,100.0d0)
3623       call reada(controlcard,'BOXY',boxysize,100.0d0)
3624       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3625       call readi(controlcard,'TUBEMOD',tubemode,0)
3626       write (iout,*) TUBEmode,"TUBEMODE"
3627       if (TUBEmode.gt.0) then
3628        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3629        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3630        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3631        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3632        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3633        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3634        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3635        buftubebot=bordtubebot+tubebufthick
3636        buftubetop=bordtubetop-tubebufthick
3637       endif
3638
3639 ! CUTOFFF ON ELECTROSTATICS
3640       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3641       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3642       write(iout,*) "R_CUT_ELE=",r_cut_ele
3643 ! Lipidic parameters
3644       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3645       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3646       if (lipthick.gt.0.0d0) then
3647        bordliptop=(boxzsize+lipthick)/2.0
3648        bordlipbot=bordliptop-lipthick
3649       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3650       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3651       buflipbot=bordlipbot+lipbufthick
3652       bufliptop=bordliptop-lipbufthick
3653       if ((lipbufthick*2.0d0).gt.lipthick) &
3654        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3655       endif !lipthick.gt.0
3656       write(iout,*) "bordliptop=",bordliptop
3657       write(iout,*) "bordlipbot=",bordlipbot
3658       write(iout,*) "bufliptop=",bufliptop
3659       write(iout,*) "buflipbot=",buflipbot
3660       write (iout,*) "SHIELD MODE",shield_mode
3661
3662 !C-------------------------
3663       minim=(index(controlcard,'MINIMIZE').gt.0)
3664       dccart=(index(controlcard,'CART').gt.0)
3665       overlapsc=(index(controlcard,'OVERLAP').gt.0)
3666       overlapsc=.not.overlapsc
3667       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3668       searchsc=.not.searchsc
3669       sideadd=(index(controlcard,'SIDEADD').gt.0)
3670       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3671       outpdb=(index(controlcard,'PDBOUT').gt.0)
3672       outmol2=(index(controlcard,'MOL2OUT').gt.0)
3673       pdbref=(index(controlcard,'PDBREF').gt.0)
3674       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3675       indpdb=index(controlcard,'PDBSTART')
3676       extconf=(index(controlcard,'EXTCONF').gt.0)
3677       call readi(controlcard,'IPRINT',iprint,0)
3678       call readi(controlcard,'MAXGEN',maxgen,10000)
3679       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3680       call readi(controlcard,"KDIAG",kdiag,0)
3681       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3682       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3683        write (iout,*) "RESCALE_MODE",rescale_mode
3684       split_ene=index(controlcard,'SPLIT_ENE').gt.0
3685       if (index(controlcard,'REGULAR').gt.0.0D0) then
3686         call reada(controlcard,'WEIDIS',weidis,0.1D0)
3687         modecalc=1
3688         refstr=.true.
3689       endif
3690       if (index(controlcard,'CHECKGRAD').gt.0) then
3691         modecalc=5
3692         if (index(controlcard,'CART').gt.0) then
3693           icheckgrad=1
3694         elseif (index(controlcard,'CARINT').gt.0) then
3695           icheckgrad=2
3696         else
3697           icheckgrad=3
3698         endif
3699       elseif (index(controlcard,'THREAD').gt.0) then
3700         modecalc=2
3701         call readi(controlcard,'THREAD',nthread,0)
3702         if (nthread.gt.0) then
3703           call reada(controlcard,'WEIDIS',weidis,0.1D0)
3704         else
3705           if (fg_rank.eq.0) &
3706           write (iout,'(a)')'A number has to follow the THREAD keyword.'
3707           stop 'Error termination in Read_Control.'
3708         endif
3709       else if (index(controlcard,'MCMA').gt.0) then
3710         modecalc=3
3711       else if (index(controlcard,'MCEE').gt.0) then
3712         modecalc=6
3713       else if (index(controlcard,'MULTCONF').gt.0) then
3714         modecalc=4
3715       else if (index(controlcard,'MAP').gt.0) then
3716         modecalc=7
3717         call readi(controlcard,'MAP',nmap,0)
3718       else if (index(controlcard,'CSA').gt.0) then
3719         modecalc=8
3720 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3721 !rc
3722 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3723 !rc
3724 !rc        modecalc=9
3725 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3726 !fmc        modecalc=10
3727       else if (index(controlcard,'SOFTREG').gt.0) then
3728         modecalc=11
3729       else if (index(controlcard,'CHECK_BOND').gt.0) then
3730         modecalc=-1
3731       else if (index(controlcard,'TEST').gt.0) then
3732         modecalc=-2
3733       else if (index(controlcard,'MD').gt.0) then
3734         modecalc=12
3735       else if (index(controlcard,'RE ').gt.0) then
3736         modecalc=14
3737       endif
3738
3739       lmuca=index(controlcard,'MUCA').gt.0
3740       call readi(controlcard,'MUCADYN',mucadyn,0)
3741       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3742       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3743        then
3744        write (iout,*) 'MUCADYN=',mucadyn
3745        write (iout,*) 'MUCASMOOTH=',muca_smooth
3746       endif
3747
3748       iscode=index(controlcard,'ONE_LETTER')
3749       indphi=index(controlcard,'PHI')
3750       indback=index(controlcard,'BACK')
3751       iranconf=index(controlcard,'RAND_CONF')
3752       i2ndstr=index(controlcard,'USE_SEC_PRED')
3753       gradout=index(controlcard,'GRADOUT').gt.0
3754       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3755       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3756       if (me.eq.king .or. .not.out1file ) &
3757         write (iout,*) "DISTCHAINMAX",distchainmax
3758
3759       if(me.eq.king.or..not.out1file) &
3760        write (iout,'(2a)') diagmeth(kdiag),&
3761         ' routine used to diagonalize matrices.'
3762       if (shield_mode.gt.0) then
3763       pi=3.141592d0
3764 !C VSolvSphere the volume of solving sphere
3765 !C      print *,pi,"pi"
3766 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3767 !C there will be no distinction between proline peptide group and normal peptide
3768 !C group in case of shielding parameters
3769       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3770       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3771       write (iout,*) VSolvSphere,VSolvSphere_div
3772 !C long axis of side chain
3773       do i=1,ntyp
3774       long_r_sidechain(i)=vbldsc0(1,i)
3775       short_r_sidechain(i)=sigma0(i)
3776       write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3777          sigma0(i)
3778       enddo
3779       buff_shield=1.0d0
3780       endif
3781       return
3782       end subroutine read_control
3783 !-----------------------------------------------------------------------------
3784       subroutine read_REMDpar
3785 !
3786 ! Read REMD settings
3787 !
3788 !       use control
3789 !       use energy
3790 !       use geometry
3791       use REMD_data
3792       use MPI_data
3793       use control_data, only:out1file
3794 !      implicit real*8 (a-h,o-z)
3795 !      include 'DIMENSIONS'
3796 !      include 'COMMON.IOUNITS'
3797 !      include 'COMMON.TIME1'
3798 !      include 'COMMON.MD'
3799       use MD_data
3800 !el #ifndef LANG0
3801 !el      include 'COMMON.LANGEVIN'
3802 !el #else
3803 !el      include 'COMMON.LANGEVIN.lang0'
3804 !el #endif
3805 !      include 'COMMON.INTERACT'
3806 !      include 'COMMON.NAMES'
3807 !      include 'COMMON.GEO'
3808 !      include 'COMMON.REMD'
3809 !      include 'COMMON.CONTROL'
3810 !      include 'COMMON.SETUP'
3811 !      character(len=80) :: ucase
3812       character(len=320) :: controlcard
3813       character(len=3200) :: controlcard1
3814       integer :: iremd_m_total
3815 !el local variables
3816       integer :: i
3817 !     real(kind=8) :: var,ene
3818
3819       if(me.eq.king.or..not.out1file) &
3820        write (iout,*) "REMD setup"
3821
3822       call card_concat(controlcard,.true.)
3823       call readi(controlcard,"NREP",nrep,3)
3824       call readi(controlcard,"NSTEX",nstex,1000)
3825       call reada(controlcard,"RETMIN",retmin,10.0d0)
3826       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3827       mremdsync=(index(controlcard,'SYNC').gt.0)
3828       call readi(controlcard,"NSYN",i_sync_step,100)
3829       restart1file=(index(controlcard,'REST1FILE').gt.0)
3830       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3831       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3832       if(max_cache_traj_use.gt.max_cache_traj) &
3833                 max_cache_traj_use=max_cache_traj
3834       if(me.eq.king.or..not.out1file) then
3835 !d       if (traj1file) then
3836 !rc caching is in testing - NTWX is not ignored
3837 !d        write (iout,*) "NTWX value is ignored"
3838 !d        write (iout,*) "  trajectory is stored to one file by master"
3839 !d        write (iout,*) "  before exchange at NSTEX intervals"
3840 !d       endif
3841        write (iout,*) "NREP= ",nrep
3842        write (iout,*) "NSTEX= ",nstex
3843        write (iout,*) "SYNC= ",mremdsync
3844        write (iout,*) "NSYN= ",i_sync_step
3845        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3846       endif
3847       remd_tlist=.false.
3848       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3849       if (index(controlcard,'TLIST').gt.0) then
3850          remd_tlist=.true.
3851          call card_concat(controlcard1,.true.)
3852          read(controlcard1,*) (remd_t(i),i=1,nrep)
3853          if(me.eq.king.or..not.out1file) &
3854           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3855       endif
3856       remd_mlist=.false.
3857       if (index(controlcard,'MLIST').gt.0) then
3858          remd_mlist=.true.
3859          call card_concat(controlcard1,.true.)
3860          read(controlcard1,*) (remd_m(i),i=1,nrep)
3861          if(me.eq.king.or..not.out1file) then
3862           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3863           iremd_m_total=0
3864           do i=1,nrep
3865            iremd_m_total=iremd_m_total+remd_m(i)
3866           enddo
3867           write (iout,*) 'Total number of replicas ',iremd_m_total
3868          endif
3869       endif
3870       if(me.eq.king.or..not.out1file) &
3871        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3872       return
3873       end subroutine read_REMDpar
3874 !-----------------------------------------------------------------------------
3875       subroutine read_MDpar
3876 !
3877 ! Read MD settings
3878 !
3879       use control_data, only: r_cut,rlamb,out1file
3880       use energy_data
3881       use geometry_data, only: pi
3882       use MPI_data
3883 !      implicit real*8 (a-h,o-z)
3884 !      include 'DIMENSIONS'
3885 !      include 'COMMON.IOUNITS'
3886 !      include 'COMMON.TIME1'
3887 !      include 'COMMON.MD'
3888       use MD_data
3889 !el #ifndef LANG0
3890 !el      include 'COMMON.LANGEVIN'
3891 !el #else
3892 !el      include 'COMMON.LANGEVIN.lang0'
3893 !el #endif
3894 !      include 'COMMON.INTERACT'
3895 !      include 'COMMON.NAMES'
3896 !      include 'COMMON.GEO'
3897 !      include 'COMMON.SETUP'
3898 !      include 'COMMON.CONTROL'
3899 !      include 'COMMON.SPLITELE'
3900 !      character(len=80) :: ucase
3901       character(len=320) :: controlcard
3902 !el local variables
3903       integer :: i,j
3904       real(kind=8) :: eta
3905
3906       call card_concat(controlcard,.true.)
3907       call readi(controlcard,"NSTEP",n_timestep,1000000)
3908       call readi(controlcard,"NTWE",ntwe,100)
3909       call readi(controlcard,"NTWX",ntwx,1000)
3910       call reada(controlcard,"DT",d_time,1.0d-1)
3911       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3912       call reada(controlcard,"DAMAX",damax,1.0d1)
3913       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3914       call readi(controlcard,"LANG",lang,0)
3915       RESPA = index(controlcard,"RESPA") .gt. 0
3916       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3917       ntime_split0=ntime_split
3918       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3919       ntime_split0=ntime_split
3920       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3921       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3922       rest = index(controlcard,"REST").gt.0
3923       tbf = index(controlcard,"TBF").gt.0
3924       usampl = index(controlcard,"USAMPL").gt.0
3925       mdpdb = index(controlcard,"MDPDB").gt.0
3926       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3927       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3928       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3929       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3930       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3931       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3932       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3933       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3934       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3935       large = index(controlcard,"LARGE").gt.0
3936       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3937       rattle = index(controlcard,"RATTLE").gt.0
3938 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3939       nset=0
3940       if(usampl) then
3941         call read_fragments
3942       endif
3943
3944       if(me.eq.king.or..not.out1file) then
3945        write (iout,*)
3946        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3947        write (iout,*)
3948        write (iout,'(a)') "The units are:"
3949        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3950        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3951         " acceleration: angstrom/(48.9 fs)**2"
3952        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3953        write (iout,*)
3954        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3955        write (iout,'(a60,f10.5,a)') &
3956         "Initial time step of numerical integration:",d_time,&
3957         " natural units"
3958        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3959        if (RESPA) then
3960         write (iout,'(2a,i4,a)') &
3961           "A-MTS algorithm used; initial time step for fast-varying",&
3962           " short-range forces split into",ntime_split," steps."
3963         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3964          r_cut," lambda",rlamb
3965        endif
3966        write (iout,'(2a,f10.5)') &
3967         "Maximum acceleration threshold to reduce the time step",&
3968         "/increase split number:",damax
3969        write (iout,'(2a,f10.5)') &
3970         "Maximum predicted energy drift to reduce the timestep",&
3971         "/increase split number:",edriftmax
3972        write (iout,'(a60,f10.5)') &
3973        "Maximum velocity threshold to reduce velocities:",dvmax
3974        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3975        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3976        if (rattle) write (iout,'(a60)') &
3977         "Rattle algorithm used to constrain the virtual bonds"
3978       endif
3979       reset_fricmat=1000
3980       if (lang.gt.0) then
3981         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3982         call reada(controlcard,"RWAT",rwat,1.4d0)
3983         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3984         surfarea=index(controlcard,"SURFAREA").gt.0
3985         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3986         if(me.eq.king.or..not.out1file)then
3987          write (iout,'(/a,$)') "Langevin dynamics calculation"
3988          if (lang.eq.1) then
3989           write (iout,'(a/)') &
3990             " with direct integration of Langevin equations"
3991          else if (lang.eq.2) then
3992           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3993          else if (lang.eq.3) then
3994           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3995          else if (lang.eq.4) then
3996           write (iout,'(a/)') " in overdamped mode"
3997          else
3998           write (iout,'(//a,i5)') &
3999             "=========== ERROR: Unknown Langevin dynamics mode:",lang
4000           stop
4001          endif
4002          write (iout,'(a60,f10.5)') "Temperature:",t_bath
4003          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4004          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4005          write (iout,'(a60,f10.5)') &
4006          "Scaling factor of the friction forces:",scal_fric
4007          if (surfarea) write (iout,'(2a,i10,a)') &
4008            "Friction coefficients will be scaled by solvent-accessible",&
4009            " surface area every",reset_fricmat," steps."
4010         endif
4011 ! Calculate friction coefficients and bounds of stochastic forces
4012         eta=6*pi*cPoise*etawat
4013         if(me.eq.king.or..not.out1file) &
4014          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4015           eta
4016 !        allocate(gamp
4017         do j=1,5 !types of molecules
4018         gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4019         stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4020         enddo
4021         allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4022         do j=1,5 !types of molecules
4023         do i=1,ntyp
4024           gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4025           stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4026         enddo
4027         enddo
4028
4029         if(me.eq.king.or..not.out1file)then
4030          write (iout,'(/2a/)') &
4031          "Radii of site types and friction coefficients and std's of",&
4032          " stochastic forces of fully exposed sites"
4033          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4034          do i=1,ntyp
4035           write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4036            gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4037          enddo
4038         endif
4039       else if (tbf) then
4040         if(me.eq.king.or..not.out1file)then
4041          write (iout,'(a)') "Berendsen bath calculation"
4042          write (iout,'(a60,f10.5)') "Temperature:",t_bath
4043          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4044          if (reset_moment) &
4045          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4046          count_reset_moment," steps"
4047          if (reset_vel) &
4048           write (iout,'(a,i10,a)') &
4049           "Velocities will be reset at random every",count_reset_vel,&
4050          " steps"
4051         endif
4052       else
4053         if(me.eq.king.or..not.out1file) &
4054          write (iout,'(a31)') "Microcanonical mode calculation"
4055       endif
4056       if(me.eq.king.or..not.out1file)then
4057        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4058        if (usampl) then
4059           write(iout,*) "MD running with constraints."
4060           write(iout,*) "Equilibration time ", eq_time, " mtus."
4061           write(iout,*) "Constraining ", nfrag," fragments."
4062           write(iout,*) "Length of each fragment, weight and q0:"
4063           do iset=1,nset
4064            write (iout,*) "Set of restraints #",iset
4065            do i=1,nfrag
4066               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4067                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4068            enddo
4069            write(iout,*) "constraints between ", npair, "fragments."
4070            write(iout,*) "constraint pairs, weights and q0:"
4071            do i=1,npair
4072             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4073                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4074            enddo
4075            write(iout,*) "angle constraints within ", nfrag_back,&
4076             "backbone fragments."
4077            write(iout,*) "fragment, weights:"
4078            do i=1,nfrag_back
4079             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4080                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4081                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4082            enddo
4083           enddo
4084         iset=mod(kolor,nset)+1
4085        endif
4086       endif
4087       if(me.eq.king.or..not.out1file) &
4088        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4089       return
4090       end subroutine read_MDpar
4091 !-----------------------------------------------------------------------------
4092       subroutine map_read
4093
4094       use map_data
4095 !      implicit real*8 (a-h,o-z)
4096 !      include 'DIMENSIONS'
4097 !      include 'COMMON.MAP'
4098 !      include 'COMMON.IOUNITS'
4099       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4100       character(len=80) :: mapcard      !,ucase
4101 !el local variables
4102       integer :: imap
4103 !     real(kind=8) :: var,ene
4104
4105       do imap=1,nmap
4106         read (inp,'(a)') mapcard
4107         mapcard=ucase(mapcard)
4108         if (index(mapcard,'PHI').gt.0) then
4109           kang(imap)=1
4110         else if (index(mapcard,'THE').gt.0) then
4111           kang(imap)=2
4112         else if (index(mapcard,'ALP').gt.0) then
4113           kang(imap)=3
4114         else if (index(mapcard,'OME').gt.0) then
4115           kang(imap)=4
4116         else
4117           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4118           stop 'Error - illegal variable spec in MAP card.'
4119         endif
4120         call readi (mapcard,'RES1',res1(imap),0)
4121         call readi (mapcard,'RES2',res2(imap),0)
4122         if (res1(imap).eq.0) then
4123           res1(imap)=res2(imap)
4124         else if (res2(imap).eq.0) then
4125           res2(imap)=res1(imap)
4126         endif
4127         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4128           write (iout,'(a)') &
4129           'Error - illegal definition of variable group in MAP.'
4130           stop 'Error - illegal definition of variable group in MAP.'
4131         endif
4132         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4133         call reada(mapcard,'TO',ang_to(imap),0.0D0)
4134         call readi(mapcard,'NSTEP',nstep(imap),0)
4135         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4136           write (iout,'(a)') &
4137            'Illegal boundary and/or step size specification in MAP.'
4138           stop 'Illegal boundary and/or step size specification in MAP.'
4139         endif
4140       enddo ! imap
4141       return
4142       end subroutine map_read
4143 !-----------------------------------------------------------------------------
4144       subroutine csaread
4145
4146       use control_data, only: vdisulf
4147       use csa_data
4148 !      implicit real*8 (a-h,o-z)
4149 !      include 'DIMENSIONS'
4150 !      include 'COMMON.IOUNITS'
4151 !      include 'COMMON.GEO'
4152 !      include 'COMMON.CSA'
4153 !      include 'COMMON.BANK'
4154 !      include 'COMMON.CONTROL'
4155 !      character(len=80) :: ucase
4156       character(len=620) :: mcmcard
4157 !el local variables
4158 !     integer :: ntf,ik,iw_pdb
4159 !     real(kind=8) :: var,ene
4160
4161       call card_concat(mcmcard,.true.)
4162
4163       call readi(mcmcard,'NCONF',nconf,50)
4164       call readi(mcmcard,'NADD',nadd,0)
4165       call readi(mcmcard,'JSTART',jstart,1)
4166       call readi(mcmcard,'JEND',jend,1)
4167       call readi(mcmcard,'NSTMAX',nstmax,500000)
4168       call readi(mcmcard,'N0',n0,1)
4169       call readi(mcmcard,'N1',n1,6)
4170       call readi(mcmcard,'N2',n2,4)
4171       call readi(mcmcard,'N3',n3,0)
4172       call readi(mcmcard,'N4',n4,0)
4173       call readi(mcmcard,'N5',n5,0)
4174       call readi(mcmcard,'N6',n6,10)
4175       call readi(mcmcard,'N7',n7,0)
4176       call readi(mcmcard,'N8',n8,0)
4177       call readi(mcmcard,'N9',n9,0)
4178       call readi(mcmcard,'N14',n14,0)
4179       call readi(mcmcard,'N15',n15,0)
4180       call readi(mcmcard,'N16',n16,0)
4181       call readi(mcmcard,'N17',n17,0)
4182       call readi(mcmcard,'N18',n18,0)
4183
4184       vdisulf=(index(mcmcard,'DYNSS').gt.0)
4185
4186       call readi(mcmcard,'NDIFF',ndiff,2)
4187       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4188       call readi(mcmcard,'IS1',is1,1)
4189       call readi(mcmcard,'IS2',is2,8)
4190       call readi(mcmcard,'NRAN0',nran0,4)
4191       call readi(mcmcard,'NRAN1',nran1,2)
4192       call readi(mcmcard,'IRR',irr,1)
4193       call readi(mcmcard,'NSEED',nseed,20)
4194       call readi(mcmcard,'NTOTAL',ntotal,10000)
4195       call reada(mcmcard,'CUT1',cut1,2.0d0)
4196       call reada(mcmcard,'CUT2',cut2,5.0d0)
4197       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4198       call readi(mcmcard,'ICMAX',icmax,3)
4199       call readi(mcmcard,'IRESTART',irestart,0)
4200 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
4201       ntbankm=0
4202 !!bankt
4203       call reada(mcmcard,'DELE',dele,20.0d0)
4204       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4205       call readi(mcmcard,'IREF',iref,0)
4206       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4207       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4208       call readi(mcmcard,'NCONF_IN',nconf_in,0)
4209       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4210       write (iout,*) "NCONF_IN",nconf_in
4211       return
4212       end subroutine csaread
4213 !-----------------------------------------------------------------------------
4214       subroutine mcmread
4215
4216       use mcm_data
4217       use control_data, only: MaxMoveType
4218       use MD_data
4219       use minim_data
4220 !      implicit real*8 (a-h,o-z)
4221 !      include 'DIMENSIONS'
4222 !      include 'COMMON.MCM'
4223 !      include 'COMMON.MCE'
4224 !      include 'COMMON.IOUNITS'
4225 !      character(len=80) :: ucase
4226       character(len=320) :: mcmcard
4227 !el local variables
4228       integer :: i
4229 !     real(kind=8) :: var,ene
4230
4231       call card_concat(mcmcard,.true.)
4232       call readi(mcmcard,'MAXACC',maxacc,100)
4233       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4234       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4235       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4236       call readi(mcmcard,'MAXREPM',maxrepm,200)
4237       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4238       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4239       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4240       call reada(mcmcard,'E_UP',e_up,5.0D0)
4241       call reada(mcmcard,'DELTE',delte,0.1D0)
4242       call readi(mcmcard,'NSWEEP',nsweep,5)
4243       call readi(mcmcard,'NSTEPH',nsteph,0)
4244       call readi(mcmcard,'NSTEPC',nstepc,0)
4245       call reada(mcmcard,'TMIN',tmin,298.0D0)
4246       call reada(mcmcard,'TMAX',tmax,298.0D0)
4247       call readi(mcmcard,'NWINDOW',nwindow,0)
4248       call readi(mcmcard,'PRINT_MC',print_mc,0)
4249       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4250       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4251       ent_read=(index(mcmcard,'ENT_READ').gt.0)
4252       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4253       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4254       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4255       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4256       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4257       if (nwindow.gt.0) then
4258         allocate(winstart(nwindow))     !!el (maxres)
4259         allocate(winend(nwindow))       !!el
4260         allocate(winlen(nwindow))       !!el
4261         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4262         do i=1,nwindow
4263           winlen(i)=winend(i)-winstart(i)+1
4264         enddo
4265       endif
4266       if (tmax.lt.tmin) tmax=tmin
4267       if (tmax.eq.tmin) then
4268         nstepc=0
4269         nsteph=0
4270       endif
4271       if (nstepc.gt.0 .and. nsteph.gt.0) then
4272         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4273         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4274       endif
4275       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4276 ! Probabilities of different move types
4277       sumpro_type(0)=0.0D0
4278       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4279       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4280       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4281       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
4282       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4283       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4284       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4285       do i=1,MaxMoveType
4286         print *,'i',i,' sumprotype',sumpro_type(i)
4287         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4288         print *,'i',i,' sumprotype',sumpro_type(i)
4289       enddo
4290       return
4291       end subroutine mcmread
4292 !-----------------------------------------------------------------------------
4293       subroutine read_minim
4294
4295       use minim_data
4296 !      implicit real*8 (a-h,o-z)
4297 !      include 'DIMENSIONS'
4298 !      include 'COMMON.MINIM'
4299 !      include 'COMMON.IOUNITS'
4300 !      character(len=80) :: ucase
4301       character(len=320) :: minimcard
4302 !el local variables
4303 !     integer :: ntf,ik,iw_pdb
4304 !     real(kind=8) :: var,ene
4305
4306       call card_concat(minimcard,.true.)
4307       call readi(minimcard,'MAXMIN',maxmin,2000)
4308       call readi(minimcard,'MAXFUN',maxfun,5000)
4309       call readi(minimcard,'MINMIN',minmin,maxmin)
4310       call readi(minimcard,'MINFUN',minfun,maxmin)
4311       call reada(minimcard,'TOLF',tolf,1.0D-2)
4312       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4313       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4314       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4315       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4316       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4317                'Options in energy minimization:'
4318       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4319        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4320        'MinMin:',MinMin,' MinFun:',MinFun,&
4321        ' TolF:',TolF,' RTolF:',RTolF
4322       return
4323       end subroutine read_minim
4324 !-----------------------------------------------------------------------------
4325       subroutine openunits
4326
4327       use MD_data, only: usampl
4328       use csa_data
4329       use MPI_data
4330       use control_data, only:out1file
4331       use control, only: getenv_loc
4332 !      implicit real*8 (a-h,o-z)
4333 !      include 'DIMENSIONS'
4334 #ifdef MPI
4335       include 'mpif.h'
4336       character(len=16) :: form,nodename
4337       integer :: nodelen,ierror,npos
4338 #endif
4339 !      include 'COMMON.SETUP'
4340 !      include 'COMMON.IOUNITS'
4341 !      include 'COMMON.MD'
4342 !      include 'COMMON.CONTROL'
4343       integer :: lenpre,lenpot,lentmp   !,ilen
4344 !el      external ilen
4345       character(len=3) :: out1file_text !,ucase
4346       character(len=3) :: ll
4347 !el      external ucase
4348 !el local variables
4349 !     integer :: ntf,ik,iw_pdb
4350 !     real(kind=8) :: var,ene
4351 !
4352 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4353       call getenv_loc("PREFIX",prefix)
4354       pref_orig = prefix
4355       call getenv_loc("POT",pot)
4356       call getenv_loc("DIRTMP",tmpdir)
4357       call getenv_loc("CURDIR",curdir)
4358       call getenv_loc("OUT1FILE",out1file_text)
4359 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4360       out1file_text=ucase(out1file_text)
4361       if (out1file_text(1:1).eq."Y") then
4362         out1file=.true.
4363       else
4364         out1file=fg_rank.gt.0
4365       endif
4366       lenpre=ilen(prefix)
4367       lenpot=ilen(pot)
4368       lentmp=ilen(tmpdir)
4369       if (lentmp.gt.0) then
4370           write (*,'(80(1h!))')
4371           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
4372           write (*,'(80(1h!))')
4373           write (*,*)"All output files will be on node /tmp directory."
4374 #ifdef MPI
4375         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4376         if (me.eq.king) then
4377           write (*,*) "The master node is ",nodename
4378         else if (fg_rank.eq.0) then
4379           write (*,*) "I am the CG slave node ",nodename
4380         else
4381           write (*,*) "I am the FG slave node ",nodename
4382         endif
4383 #endif
4384         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4385         lenpre = lentmp+lenpre+1
4386       endif
4387       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4388 ! Get the names and open the input files
4389 #if defined(WINIFL) || defined(WINPGI)
4390       open(1,file=pref_orig(:ilen(pref_orig))// &
4391         '.inp',status='old',readonly,shared)
4392        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4393 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4394 ! Get parameter filenames and open the parameter files.
4395       call getenv_loc('BONDPAR',bondname)
4396       open (ibond,file=bondname,status='old',readonly,shared)
4397       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4398       open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4399       call getenv_loc('THETPAR',thetname)
4400       open (ithep,file=thetname,status='old',readonly,shared)
4401       call getenv_loc('ROTPAR',rotname)
4402       open (irotam,file=rotname,status='old',readonly,shared)
4403       call getenv_loc('TORPAR',torname)
4404       open (itorp,file=torname,status='old',readonly,shared)
4405       call getenv_loc('TORDPAR',tordname)
4406       open (itordp,file=tordname,status='old',readonly,shared)
4407       call getenv_loc('FOURIER',fouriername)
4408       open (ifourier,file=fouriername,status='old',readonly,shared)
4409       call getenv_loc('ELEPAR',elename)
4410       open (ielep,file=elename,status='old',readonly,shared)
4411       call getenv_loc('SIDEPAR',sidename)
4412       open (isidep,file=sidename,status='old',readonly,shared)
4413
4414       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4415       open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4416       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4417       open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4418       call getenv_loc('TORPAR_NUCL',torname_nucl)
4419       open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4420       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4421       open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4422       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4423       open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4424
4425
4426 #elif (defined CRAY) || (defined AIX)
4427       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4428         action='read')
4429 !      print *,"Processor",myrank," opened file 1"
4430       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4431 !      print *,"Processor",myrank," opened file 9"
4432 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4433 ! Get parameter filenames and open the parameter files.
4434       call getenv_loc('BONDPAR',bondname)
4435       open (ibond,file=bondname,status='old',action='read')
4436       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4437       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4438
4439 !      print *,"Processor",myrank," opened file IBOND"
4440       call getenv_loc('THETPAR',thetname)
4441       open (ithep,file=thetname,status='old',action='read')
4442 !      print *,"Processor",myrank," opened file ITHEP"
4443       call getenv_loc('ROTPAR',rotname)
4444       open (irotam,file=rotname,status='old',action='read')
4445 !      print *,"Processor",myrank," opened file IROTAM"
4446       call getenv_loc('TORPAR',torname)
4447       open (itorp,file=torname,status='old',action='read')
4448 !      print *,"Processor",myrank," opened file ITORP"
4449       call getenv_loc('TORDPAR',tordname)
4450       open (itordp,file=tordname,status='old',action='read')
4451 !      print *,"Processor",myrank," opened file ITORDP"
4452       call getenv_loc('SCCORPAR',sccorname)
4453       open (isccor,file=sccorname,status='old',action='read')
4454 !      print *,"Processor",myrank," opened file ISCCOR"
4455       call getenv_loc('FOURIER',fouriername)
4456       open (ifourier,file=fouriername,status='old',action='read')
4457 !      print *,"Processor",myrank," opened file IFOURIER"
4458       call getenv_loc('ELEPAR',elename)
4459       open (ielep,file=elename,status='old',action='read')
4460 !      print *,"Processor",myrank," opened file IELEP"
4461       call getenv_loc('SIDEPAR',sidename)
4462       open (isidep,file=sidename,status='old',action='read')
4463
4464       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4465       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4466       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4467       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4468       call getenv_loc('TORPAR_NUCL',torname_nucl)
4469       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4470       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4471       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4472       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4473       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4474
4475       call getenv_loc('LIPTRANPAR',liptranname)
4476       open (iliptranpar,file=liptranname,status='old',action='read')
4477       call getenv_loc('TUBEPAR',tubename)
4478       open (itube,file=tubename,status='old',action='read')
4479       call getenv_loc('IONPAR',ionname)
4480       open (iion,file=ionname,status='old',action='read')
4481
4482 !      print *,"Processor",myrank," opened file ISIDEP"
4483 !      print *,"Processor",myrank," opened parameter files"
4484 #elif (defined G77)
4485       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4486       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4487 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4488 ! Get parameter filenames and open the parameter files.
4489       call getenv_loc('BONDPAR',bondname)
4490       open (ibond,file=bondname,status='old')
4491       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4492       open (ibond_nucl,file=bondname_nucl,status='old')
4493
4494       call getenv_loc('THETPAR',thetname)
4495       open (ithep,file=thetname,status='old')
4496       call getenv_loc('ROTPAR',rotname)
4497       open (irotam,file=rotname,status='old')
4498       call getenv_loc('TORPAR',torname)
4499       open (itorp,file=torname,status='old')
4500       call getenv_loc('TORDPAR',tordname)
4501       open (itordp,file=tordname,status='old')
4502       call getenv_loc('SCCORPAR',sccorname)
4503       open (isccor,file=sccorname,status='old')
4504       call getenv_loc('FOURIER',fouriername)
4505       open (ifourier,file=fouriername,status='old')
4506       call getenv_loc('ELEPAR',elename)
4507       open (ielep,file=elename,status='old')
4508       call getenv_loc('SIDEPAR',sidename)
4509       open (isidep,file=sidename,status='old')
4510
4511       open (ithep_nucl,file=thetname_nucl,status='old')
4512       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4513       open (irotam_nucl,file=rotname_nucl,status='old')
4514       call getenv_loc('TORPAR_NUCL',torname_nucl)
4515       open (itorp_nucl,file=torname_nucl,status='old')
4516       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4517       open (itordp_nucl,file=tordname_nucl,status='old')
4518       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4519       open (isidep_nucl,file=sidename_nucl,status='old')
4520
4521       call getenv_loc('LIPTRANPAR',liptranname)
4522       open (iliptranpar,file=liptranname,status='old')
4523       call getenv_loc('TUBEPAR',tubename)
4524       open (itube,file=tubename,status='old')
4525       call getenv_loc('IONPAR',ionname)
4526       open (iion,file=ionname,status='old')
4527 #else
4528       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4529         readonly)
4530        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4531 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4532 ! Get parameter filenames and open the parameter files.
4533       call getenv_loc('BONDPAR',bondname)
4534       open (ibond,file=bondname,status='old',action='read')
4535       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4536       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4537       call getenv_loc('THETPAR',thetname)
4538       open (ithep,file=thetname,status='old',action='read')
4539       call getenv_loc('ROTPAR',rotname)
4540       open (irotam,file=rotname,status='old',action='read')
4541       call getenv_loc('TORPAR',torname)
4542       open (itorp,file=torname,status='old',action='read')
4543       call getenv_loc('TORDPAR',tordname)
4544       open (itordp,file=tordname,status='old',action='read')
4545       call getenv_loc('SCCORPAR',sccorname)
4546       open (isccor,file=sccorname,status='old',action='read')
4547 #ifndef CRYST_THETA
4548       call getenv_loc('THETPARPDB',thetname_pdb)
4549       print *,"thetname_pdb ",thetname_pdb
4550       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4551       print *,ithep_pdb," opened"
4552 #endif
4553       call getenv_loc('FOURIER',fouriername)
4554       open (ifourier,file=fouriername,status='old',readonly)
4555       call getenv_loc('ELEPAR',elename)
4556       open (ielep,file=elename,status='old',readonly)
4557       call getenv_loc('SIDEPAR',sidename)
4558       open (isidep,file=sidename,status='old',readonly)
4559
4560       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4561       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4562       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4563       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4564       call getenv_loc('TORPAR_NUCL',torname_nucl)
4565       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4566       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4567       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4568       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4569       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4570       call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4571       open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4572       call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4573       open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4574       call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4575       open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4576       call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4577       open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4578
4579
4580       call getenv_loc('LIPTRANPAR',liptranname)
4581       open (iliptranpar,file=liptranname,status='old',action='read')
4582       call getenv_loc('TUBEPAR',tubename)
4583       open (itube,file=tubename,status='old',action='read')
4584       call getenv_loc('IONPAR',ionname)
4585       open (iion,file=ionname,status='old',action='read')
4586
4587 #ifndef CRYST_SC
4588       call getenv_loc('ROTPARPDB',rotname_pdb)
4589       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4590 #endif
4591 #endif
4592       call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4593 #if defined(WINIFL) || defined(WINPGI)
4594       open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4595 #elif (defined CRAY)  || (defined AIX)
4596       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4597 #elif (defined G77)
4598       open (iscpp_nucl,file=scpname_nucl,status='old')
4599 #else
4600       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4601 #endif
4602
4603 #ifndef OLDSCP
4604 !
4605 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4606 ! Use -DOLDSCP to use hard-coded constants instead.
4607 !
4608       call getenv_loc('SCPPAR',scpname)
4609 #if defined(WINIFL) || defined(WINPGI)
4610       open (iscpp,file=scpname,status='old',readonly,shared)
4611 #elif (defined CRAY)  || (defined AIX)
4612       open (iscpp,file=scpname,status='old',action='read')
4613 #elif (defined G77)
4614       open (iscpp,file=scpname,status='old')
4615 #else
4616       open (iscpp,file=scpname,status='old',action='read')
4617 #endif
4618 #endif
4619       call getenv_loc('PATTERN',patname)
4620 #if defined(WINIFL) || defined(WINPGI)
4621       open (icbase,file=patname,status='old',readonly,shared)
4622 #elif (defined CRAY)  || (defined AIX)
4623       open (icbase,file=patname,status='old',action='read')
4624 #elif (defined G77)
4625       open (icbase,file=patname,status='old')
4626 #else
4627       open (icbase,file=patname,status='old',action='read')
4628 #endif
4629 #ifdef MPI
4630 ! Open output file only for CG processes
4631 !      print *,"Processor",myrank," fg_rank",fg_rank
4632       if (fg_rank.eq.0) then
4633
4634       if (nodes.eq.1) then
4635         npos=3
4636       else
4637         npos = dlog10(dfloat(nodes-1))+1
4638       endif
4639       if (npos.lt.3) npos=3
4640       write (liczba,'(i1)') npos
4641       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4642         //')'
4643       write (liczba,form) me
4644       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4645         liczba(:ilen(liczba))
4646       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4647         //'.int'
4648       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4649         //'.pdb'
4650       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4651         liczba(:ilen(liczba))//'.mol2'
4652       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4653         liczba(:ilen(liczba))//'.stat'
4654       if (lentmp.gt.0) &
4655         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4656             //liczba(:ilen(liczba))//'.stat')
4657       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4658         //'.rst'
4659       if(usampl) then
4660           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4661        liczba(:ilen(liczba))//'.const'
4662       endif
4663
4664       endif
4665 #else
4666       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4667       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4668       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4669       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4670       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4671       if (lentmp.gt.0) &
4672         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4673           '.stat')
4674       rest2name=prefix(:ilen(prefix))//'.rst'
4675       if(usampl) then
4676          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4677       endif
4678 #endif
4679 #if defined(AIX) || defined(PGI)
4680       if (me.eq.king .or. .not. out1file) &
4681          open(iout,file=outname,status='unknown')
4682 #ifdef DEBUG
4683       if (fg_rank.gt.0) then
4684         write (liczba,'(i3.3)') myrank/nfgtasks
4685         write (ll,'(bz,i3.3)') fg_rank
4686         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4687          status='unknown')
4688       endif
4689 #endif
4690       if(me.eq.king) then
4691        open(igeom,file=intname,status='unknown',position='append')
4692        open(ipdb,file=pdbname,status='unknown')
4693        open(imol2,file=mol2name,status='unknown')
4694        open(istat,file=statname,status='unknown',position='append')
4695       else
4696 !1out       open(iout,file=outname,status='unknown')
4697       endif
4698 #else
4699       if (me.eq.king .or. .not.out1file) &
4700           open(iout,file=outname,status='unknown')
4701 #ifdef DEBUG
4702       if (fg_rank.gt.0) then
4703         write (liczba,'(i3.3)') myrank/nfgtasks
4704         write (ll,'(bz,i3.3)') fg_rank
4705         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4706          status='unknown')
4707       endif
4708 #endif
4709       if(me.eq.king) then
4710        open(igeom,file=intname,status='unknown',access='append')
4711        open(ipdb,file=pdbname,status='unknown')
4712        open(imol2,file=mol2name,status='unknown')
4713        open(istat,file=statname,status='unknown',access='append')
4714       else
4715 !1out       open(iout,file=outname,status='unknown')
4716       endif
4717 #endif
4718       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4719       csa_seed=prefix(:lenpre)//'.CSA.seed'
4720       csa_history=prefix(:lenpre)//'.CSA.history'
4721       csa_bank=prefix(:lenpre)//'.CSA.bank'
4722       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4723       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4724       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4725 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4726       csa_int=prefix(:lenpre)//'.int'
4727       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4728       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4729       csa_in=prefix(:lenpre)//'.CSA.in'
4730 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4731 ! Write file names
4732       if (me.eq.king)then
4733       write (iout,'(80(1h-))')
4734       write (iout,'(30x,a)') "FILE ASSIGNMENT"
4735       write (iout,'(80(1h-))')
4736       write (iout,*) "Input file                      : ",&
4737         pref_orig(:ilen(pref_orig))//'.inp'
4738       write (iout,*) "Output file                     : ",&
4739         outname(:ilen(outname))
4740       write (iout,*)
4741       write (iout,*) "Sidechain potential file        : ",&
4742         sidename(:ilen(sidename))
4743 #ifndef OLDSCP
4744       write (iout,*) "SCp potential file              : ",&
4745         scpname(:ilen(scpname))
4746 #endif
4747       write (iout,*) "Electrostatic potential file    : ",&
4748         elename(:ilen(elename))
4749       write (iout,*) "Cumulant coefficient file       : ",&
4750         fouriername(:ilen(fouriername))
4751       write (iout,*) "Torsional parameter file        : ",&
4752         torname(:ilen(torname))
4753       write (iout,*) "Double torsional parameter file : ",&
4754         tordname(:ilen(tordname))
4755       write (iout,*) "SCCOR parameter file : ",&
4756         sccorname(:ilen(sccorname))
4757       write (iout,*) "Bond & inertia constant file    : ",&
4758         bondname(:ilen(bondname))
4759       write (iout,*) "Bending parameter file          : ",&
4760         thetname(:ilen(thetname))
4761       write (iout,*) "Rotamer parameter file          : ",&
4762         rotname(:ilen(rotname))
4763 !el----
4764 #ifndef CRYST_THETA
4765       write (iout,*) "Thetpdb parameter file          : ",&
4766         thetname_pdb(:ilen(thetname_pdb))
4767 #endif
4768 !el
4769       write (iout,*) "Threading database              : ",&
4770         patname(:ilen(patname))
4771       if (lentmp.ne.0) &
4772       write (iout,*)" DIRTMP                          : ",&
4773         tmpdir(:lentmp)
4774       write (iout,'(80(1h-))')
4775       endif
4776       return
4777       end subroutine openunits
4778 !-----------------------------------------------------------------------------
4779       subroutine readrst
4780
4781       use geometry_data, only: nres,dc
4782       use MD_data
4783 !      implicit real*8 (a-h,o-z)
4784 !      include 'DIMENSIONS'
4785 !      include 'COMMON.CHAIN'
4786 !      include 'COMMON.IOUNITS'
4787 !      include 'COMMON.MD'
4788 !el local variables
4789       integer ::i,j
4790 !     real(kind=8) :: var,ene
4791
4792       open(irest2,file=rest2name,status='unknown')
4793       read(irest2,*) totT,EK,potE,totE,t_bath
4794       totTafm=totT
4795 !      do i=1,2*nres
4796 ! AL 4/17/17: Now reading d_t(0,:) too
4797       do i=0,2*nres
4798          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4799       enddo
4800 !      do i=1,2*nres
4801 ! AL 4/17/17: Now reading d_c(0,:) too
4802       do i=0,2*nres
4803          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4804       enddo
4805       if(usampl) then
4806              read (irest2,*) iset
4807       endif
4808       close(irest2)
4809       return
4810       end subroutine readrst
4811 !-----------------------------------------------------------------------------
4812       subroutine read_fragments
4813
4814       use energy_data
4815 !      use geometry
4816       use control_data, only:out1file
4817       use MD_data
4818       use MPI_data
4819 !      implicit real*8 (a-h,o-z)
4820 !      include 'DIMENSIONS'
4821 #ifdef MPI
4822       include 'mpif.h'
4823 #endif
4824 !      include 'COMMON.SETUP'
4825 !      include 'COMMON.CHAIN'
4826 !      include 'COMMON.IOUNITS'
4827 !      include 'COMMON.MD'
4828 !      include 'COMMON.CONTROL'
4829 !el local variables
4830       integer :: i
4831 !     real(kind=8) :: var,ene
4832
4833       read(inp,*) nset,nfrag,npair,nfrag_back
4834
4835 !el from module energy
4836 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
4837       if(.not.allocated(wfrag_back)) then
4838         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4839         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4840
4841         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4842         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
4843
4844         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4845         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
4846       endif
4847
4848       if(me.eq.king.or..not.out1file) &
4849        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4850         " nfrag_back",nfrag_back
4851       do iset=1,nset
4852          read(inp,*) mset(iset)
4853        do i=1,nfrag
4854          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4855            qinfrag(i,iset)
4856          if(me.eq.king.or..not.out1file) &
4857           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4858            ifrag(2,i,iset), qinfrag(i,iset)
4859        enddo
4860        do i=1,npair
4861         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4862           qinpair(i,iset)
4863         if(me.eq.king.or..not.out1file) &
4864          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4865           ipair(2,i,iset), qinpair(i,iset)
4866        enddo
4867        do i=1,nfrag_back
4868         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4869            wfrag_back(3,i,iset),&
4870            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4871         if(me.eq.king.or..not.out1file) &
4872          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4873          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4874        enddo
4875       enddo
4876       return
4877       end subroutine read_fragments
4878 !-----------------------------------------------------------------------------
4879 ! shift.F       io_csa
4880 !-----------------------------------------------------------------------------
4881       subroutine csa_read
4882
4883       use csa_data
4884 !      implicit real*8 (a-h,o-z)
4885 !      include 'DIMENSIONS'
4886 !      include 'COMMON.CSA'
4887 !      include 'COMMON.BANK'
4888 !      include 'COMMON.IOUNITS'
4889 !el local variables
4890 !     integer :: ntf,ik,iw_pdb
4891 !     real(kind=8) :: var,ene
4892
4893       open(icsa_in,file=csa_in,status="old",err=100)
4894        read(icsa_in,*) nconf
4895        read(icsa_in,*) jstart,jend
4896        read(icsa_in,*) nstmax
4897        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4898        read(icsa_in,*) nran0,nran1,irr
4899        read(icsa_in,*) nseed
4900        read(icsa_in,*) ntotal,cut1,cut2
4901        read(icsa_in,*) estop
4902        read(icsa_in,*) icmax,irestart
4903        read(icsa_in,*) ntbankm,dele,difcut
4904        read(icsa_in,*) iref,rmscut,pnccut
4905        read(icsa_in,*) ndiff
4906       close(icsa_in)
4907
4908       return
4909
4910  100  continue
4911       return
4912       end subroutine csa_read
4913 !-----------------------------------------------------------------------------
4914       subroutine initial_write
4915
4916       use csa_data
4917 !      implicit real*8 (a-h,o-z)
4918 !      include 'DIMENSIONS'
4919 !      include 'COMMON.CSA'
4920 !      include 'COMMON.BANK'
4921 !      include 'COMMON.IOUNITS'
4922 !el local variables
4923 !     integer :: ntf,ik,iw_pdb
4924 !     real(kind=8) :: var,ene
4925
4926       open(icsa_seed,file=csa_seed,status="unknown")
4927        write(icsa_seed,*) "seed"
4928       close(31)
4929 #if defined(AIX) || defined(PGI)
4930        open(icsa_history,file=csa_history,status="unknown",&
4931         position="append")
4932 #else
4933        open(icsa_history,file=csa_history,status="unknown",&
4934         access="append")
4935 #endif
4936        write(icsa_history,*) nconf
4937        write(icsa_history,*) jstart,jend
4938        write(icsa_history,*) nstmax
4939        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4940        write(icsa_history,*) nran0,nran1,irr
4941        write(icsa_history,*) nseed
4942        write(icsa_history,*) ntotal,cut1,cut2
4943        write(icsa_history,*) estop
4944        write(icsa_history,*) icmax,irestart
4945        write(icsa_history,*) ntbankm,dele,difcut
4946        write(icsa_history,*) iref,rmscut,pnccut
4947        write(icsa_history,*) ndiff
4948
4949        write(icsa_history,*)
4950       close(icsa_history)
4951
4952       open(icsa_bank1,file=csa_bank1,status="unknown")
4953        write(icsa_bank1,*) 0
4954       close(icsa_bank1)
4955
4956       return
4957       end subroutine initial_write
4958 !-----------------------------------------------------------------------------
4959       subroutine restart_write
4960
4961       use csa_data
4962 !      implicit real*8 (a-h,o-z)
4963 !      include 'DIMENSIONS'
4964 !      include 'COMMON.IOUNITS'
4965 !      include 'COMMON.CSA'
4966 !      include 'COMMON.BANK'
4967 !el local variables
4968 !     integer :: ntf,ik,iw_pdb
4969 !     real(kind=8) :: var,ene
4970
4971 #if defined(AIX) || defined(PGI)
4972        open(icsa_history,file=csa_history,position="append")
4973 #else
4974        open(icsa_history,file=csa_history,access="append")
4975 #endif
4976        write(icsa_history,*)
4977        write(icsa_history,*) "This is restart"
4978        write(icsa_history,*)
4979        write(icsa_history,*) nconf
4980        write(icsa_history,*) jstart,jend
4981        write(icsa_history,*) nstmax
4982        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4983        write(icsa_history,*) nran0,nran1,irr
4984        write(icsa_history,*) nseed
4985        write(icsa_history,*) ntotal,cut1,cut2
4986        write(icsa_history,*) estop
4987        write(icsa_history,*) icmax,irestart
4988        write(icsa_history,*) ntbankm,dele,difcut
4989        write(icsa_history,*) iref,rmscut,pnccut
4990        write(icsa_history,*) ndiff
4991        write(icsa_history,*)
4992        write(icsa_history,*) "irestart is: ", irestart
4993
4994        write(icsa_history,*)
4995       close(icsa_history)
4996
4997       return
4998       end subroutine restart_write
4999 !-----------------------------------------------------------------------------
5000 ! test.F
5001 !-----------------------------------------------------------------------------
5002       subroutine write_pdb(npdb,titelloc,ee)
5003
5004 !      implicit real*8 (a-h,o-z)
5005 !      include 'DIMENSIONS'
5006 !      include 'COMMON.IOUNITS'
5007       character(len=50) :: titelloc1
5008       character*(*) :: titelloc
5009       character(len=3) :: zahl
5010       character(len=5) :: liczba5
5011       real(kind=8) :: ee
5012       integer :: npdb   !,ilen
5013 !el      external ilen
5014 !el local variables
5015       integer :: lenpre
5016 !     real(kind=8) :: var,ene
5017
5018       titelloc1=titelloc
5019       lenpre=ilen(prefix)
5020       if (npdb.lt.1000) then
5021        call numstr(npdb,zahl)
5022        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5023       else
5024         if (npdb.lt.10000) then
5025          write(liczba5,'(i1,i4)') 0,npdb
5026         else
5027          write(liczba5,'(i5)') npdb
5028         endif
5029         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5030       endif
5031       call pdbout(ee,titelloc1,ipdb)
5032       close(ipdb)
5033       return
5034       end subroutine write_pdb
5035 !-----------------------------------------------------------------------------
5036 ! thread.F
5037 !-----------------------------------------------------------------------------
5038       subroutine write_thread_summary
5039 ! Thread the sequence through a database of known structures
5040       use control_data, only: refstr
5041 !      use geometry
5042       use energy_data, only: n_ene_comp
5043       use compare_data
5044 !      implicit real*8 (a-h,o-z)
5045 !      include 'DIMENSIONS'
5046 #ifdef MPI
5047       use MPI_data      !include 'COMMON.INFO'
5048       include 'mpif.h'
5049 #endif
5050 !      include 'COMMON.CONTROL'
5051 !      include 'COMMON.CHAIN'
5052 !      include 'COMMON.DBASE'
5053 !      include 'COMMON.INTERACT'
5054 !      include 'COMMON.VAR'
5055 !      include 'COMMON.THREAD'
5056 !      include 'COMMON.FFIELD'
5057 !      include 'COMMON.SBRIDGE'
5058 !      include 'COMMON.HEADER'
5059 !      include 'COMMON.NAMES'
5060 !      include 'COMMON.IOUNITS'
5061 !      include 'COMMON.TIME1'
5062
5063       integer,dimension(maxthread) :: ip
5064       real(kind=8),dimension(0:n_ene) :: energia
5065 !el local variables
5066       integer :: i,j,ii,jj,ipj,ik,kk,ist
5067       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5068
5069       write (iout,'(30x,a/)') &
5070        '  *********** Summary threading statistics ************'
5071       write (iout,'(a)') 'Initial energies:'
5072       write (iout,'(a4,2x,a12,14a14,3a8)') &
5073        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5074        'RMSnat','NatCONT','NNCONT','RMS'
5075 ! Energy sort patterns
5076       do i=1,nthread
5077         ip(i)=i
5078       enddo
5079       do i=1,nthread-1
5080         enet=ener(n_ene-1,ip(i))
5081         jj=i
5082         do j=i+1,nthread
5083           if (ener(n_ene-1,ip(j)).lt.enet) then
5084             jj=j
5085             enet=ener(n_ene-1,ip(j))
5086           endif
5087         enddo
5088         if (jj.ne.i) then
5089           ipj=ip(jj)
5090           ip(jj)=ip(i)
5091           ip(i)=ipj
5092         endif
5093       enddo
5094       do ik=1,nthread
5095         i=ip(ik)
5096         ii=ipatt(1,i)
5097         ist=nres_base(2,ii)+ipatt(2,i)
5098         do kk=1,n_ene_comp
5099           energia(i)=ener0(kk,i)
5100         enddo
5101         etot=ener0(n_ene_comp+1,i)
5102         rmsnat=ener0(n_ene_comp+2,i)
5103         rms=ener0(n_ene_comp+3,i)
5104         frac=ener0(n_ene_comp+4,i)
5105         frac_nn=ener0(n_ene_comp+5,i)
5106
5107         if (refstr) then
5108         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5109         i,str_nam(ii),ist+1,&
5110         (energia(print_order(kk)),kk=1,nprint_ene),&
5111         etot,rmsnat,frac,frac_nn,rms
5112         else
5113         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5114         i,str_nam(ii),ist+1,&
5115         (energia(print_order(kk)),kk=1,nprint_ene),etot
5116         endif
5117       enddo
5118       write (iout,'(//a)') 'Final energies:'
5119       write (iout,'(a4,2x,a12,17a14,3a8)') &
5120        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5121        'RMSnat','NatCONT','NNCONT','RMS'
5122       do ik=1,nthread
5123         i=ip(ik)
5124         ii=ipatt(1,i)
5125         ist=nres_base(2,ii)+ipatt(2,i)
5126         do kk=1,n_ene_comp
5127           energia(kk)=ener(kk,ik)
5128         enddo
5129         etot=ener(n_ene_comp+1,i)
5130         rmsnat=ener(n_ene_comp+2,i)
5131         rms=ener(n_ene_comp+3,i)
5132         frac=ener(n_ene_comp+4,i)
5133         frac_nn=ener(n_ene_comp+5,i)
5134         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5135         i,str_nam(ii),ist+1,&
5136         (energia(print_order(kk)),kk=1,nprint_ene),&
5137         etot,rmsnat,frac,frac_nn,rms
5138       enddo
5139       write (iout,'(/a/)') 'IEXAM array:'
5140       write (iout,'(i5)') nexcl
5141       do i=1,nexcl
5142         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5143       enddo
5144       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5145        'Max. time for threading step ',max_time_for_thread,&
5146        'Average time for threading step: ',ave_time_for_thread
5147       return
5148       end subroutine write_thread_summary
5149 !-----------------------------------------------------------------------------
5150       subroutine write_stat_thread(ithread,ipattern,ist)
5151
5152       use energy_data, only: n_ene_comp
5153       use compare_data
5154 !      implicit real*8 (a-h,o-z)
5155 !      include "DIMENSIONS"
5156 !      include "COMMON.CONTROL"
5157 !      include "COMMON.IOUNITS"
5158 !      include "COMMON.THREAD"
5159 !      include "COMMON.FFIELD"
5160 !      include "COMMON.DBASE"
5161 !      include "COMMON.NAMES"
5162       real(kind=8),dimension(0:n_ene) :: energia
5163 !el local variables
5164       integer :: ithread,ipattern,ist,i
5165       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5166
5167 #if defined(AIX) || defined(PGI)
5168       open(istat,file=statname,position='append')
5169 #else
5170       open(istat,file=statname,access='append')
5171 #endif
5172       do i=1,n_ene_comp
5173         energia(i)=ener(i,ithread)
5174       enddo
5175       etot=ener(n_ene_comp+1,ithread)
5176       rmsnat=ener(n_ene_comp+2,ithread)
5177       rms=ener(n_ene_comp+3,ithread)
5178       frac=ener(n_ene_comp+4,ithread)
5179       frac_nn=ener(n_ene_comp+5,ithread)
5180       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5181         ithread,str_nam(ipattern),ist+1,&
5182         (energia(print_order(i)),i=1,nprint_ene),&
5183         etot,rmsnat,frac,frac_nn,rms
5184       close (istat)
5185       return
5186       end subroutine write_stat_thread
5187 !-----------------------------------------------------------------------------
5188 #endif
5189 !-----------------------------------------------------------------------------
5190       end module io_config