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       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2554       bad(:,:)=0.0D0
2555
2556 #ifdef OLDSCP
2557 !
2558 ! Define the SC-p interaction constants (hard-coded; old style)
2559 !
2560       do i=1,ntyp
2561 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2562 ! helix formation)
2563 !       aad(i,1)=0.3D0*4.0D0**12
2564 ! Following line for constants currently implemented
2565 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2566         aad(i,1)=1.5D0*4.0D0**12
2567 !       aad(i,1)=0.17D0*5.6D0**12
2568         aad(i,2)=aad(i,1)
2569 ! "Soft" SC-p repulsion
2570         bad(i,1)=0.0D0
2571 ! Following line for constants currently implemented
2572 !       aad(i,1)=0.3D0*4.0D0**6
2573 ! "Hard" SC-p repulsion
2574         bad(i,1)=3.0D0*4.0D0**6
2575 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2576         bad(i,2)=bad(i,1)
2577 !       aad(i,1)=0.0D0
2578 !       aad(i,2)=0.0D0
2579 !       bad(i,1)=1228.8D0
2580 !       bad(i,2)=1228.8D0
2581       enddo
2582 #else
2583 !
2584 ! 8/9/01 Read the SC-p interaction constants from file
2585 !
2586       do i=1,ntyp
2587         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2588       enddo
2589       do i=1,ntyp
2590         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2591         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2592         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2593         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2594       enddo
2595 !      lprint=.true.
2596       if (lprint) then
2597         write (iout,*) "Parameters of SC-p interactions:"
2598         do i=1,ntyp
2599           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2600            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2601         enddo
2602       endif
2603 !      lprint=.false.
2604 #endif
2605       allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2606
2607       do i=1,ntyp_molec(2)
2608         read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2609       enddo
2610       do i=1,ntyp_molec(2)
2611         aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2612         bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2613       enddo
2614       r0pp=1.12246204830937298142*5.16158
2615       epspp=4.95713/4
2616       AEES=108.661
2617       BEES=0.433246
2618
2619 !
2620 ! Define the constants of the disulfide bridge
2621 !
2622       ebr=-5.50D0
2623 !
2624 ! Old arbitrary potential - commented out.
2625 !
2626 !      dbr= 4.20D0
2627 !      fbr= 3.30D0
2628 !
2629 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2630 ! energy surface of diethyl disulfide.
2631 ! A. Liwo and U. Kozlowska, 11/24/03
2632 !
2633       D0CM = 3.78d0
2634       AKCM = 15.1d0
2635       AKTH = 11.0d0
2636       AKCT = 12.0d0
2637       V1SS =-1.08d0
2638       V2SS = 7.61d0
2639       V3SS = 13.7d0
2640 !      akcm=0.0d0
2641 !      akth=0.0d0
2642 !      akct=0.0d0
2643 !      v1ss=0.0d0
2644 !      v2ss=0.0d0
2645 !      v3ss=0.0d0
2646
2647       if(me.eq.king) then
2648       write (iout,'(/a)') "Disulfide bridge parameters:"
2649       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2650       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2651       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2652       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2653         ' v3ss:',v3ss
2654       endif
2655       return
2656   111 write (iout,*) "Error reading bending energy parameters."
2657       goto 999
2658   112 write (iout,*) "Error reading rotamer energy parameters."
2659       goto 999
2660   113 write (iout,*) "Error reading torsional energy parameters."
2661       goto 999
2662   114 write (iout,*) "Error reading double torsional energy parameters."
2663       goto 999
2664   115 write (iout,*) &
2665         "Error reading cumulant (multibody energy) parameters."
2666       goto 999
2667   116 write (iout,*) "Error reading electrostatic energy parameters."
2668       goto 999
2669   117 write (iout,*) "Error reading side chain interaction parameters."
2670       goto 999
2671   118 write (iout,*) "Error reading SCp interaction parameters."
2672       goto 999
2673   119 write (iout,*) "Error reading SCCOR parameters"
2674   999 continue
2675 #ifdef MPI
2676       call MPI_Finalize(Ierror)
2677 #endif
2678       stop
2679       return
2680       end subroutine parmread
2681 #endif
2682 !-----------------------------------------------------------------------------
2683 ! printmat.f
2684 !-----------------------------------------------------------------------------
2685       subroutine printmat(ldim,m,n,iout,key,a)
2686
2687       integer :: n,ldim
2688       character(len=3),dimension(n) :: key
2689       real(kind=8),dimension(ldim,n) :: a
2690 !el local variables
2691       integer :: i,j,k,m,iout,nlim
2692
2693       do 1 i=1,n,8
2694       nlim=min0(i+7,n)
2695       write (iout,1000) (key(k),k=i,nlim)
2696       write (iout,1020)
2697  1000 format (/5x,8(6x,a3))
2698  1020 format (/80(1h-)/)
2699       do 2 j=1,n
2700       write (iout,1010) key(j),(a(j,k),k=i,nlim)
2701     2 continue
2702     1 continue
2703  1010 format (a3,2x,8(f9.4))
2704       return
2705       end subroutine printmat
2706 !-----------------------------------------------------------------------------
2707 ! readpdb.F
2708 !-----------------------------------------------------------------------------
2709       subroutine readpdb
2710 ! Read the PDB file and convert the peptide geometry into virtual-chain
2711 ! geometry.
2712       use geometry_data
2713       use energy_data, only: itype,istype
2714       use control_data
2715       use compare_data
2716       use MPI_data
2717       use control, only: rescode,sugarcode
2718 !      implicit real*8 (a-h,o-z)
2719 !      include 'DIMENSIONS'
2720 !      include 'COMMON.LOCAL'
2721 !      include 'COMMON.VAR'
2722 !      include 'COMMON.CHAIN'
2723 !      include 'COMMON.INTERACT'
2724 !      include 'COMMON.IOUNITS'
2725 !      include 'COMMON.GEO'
2726 !      include 'COMMON.NAMES'
2727 !      include 'COMMON.CONTROL'
2728 !      include 'COMMON.DISTFIT'
2729 !      include 'COMMON.SETUP'
2730       integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2731 !        ishift_pdb
2732       logical :: lprn=.true.,fail
2733       real(kind=8),dimension(3) :: e1,e2,e3
2734       real(kind=8) :: dcj,efree_temp
2735       character(len=3) :: seq,res
2736       character(len=5) :: atom
2737       character(len=80) :: card
2738       real(kind=8),dimension(3,20) :: sccor
2739       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin       !rescode,
2740       integer :: isugar,molecprev,firstion
2741       character*1 :: sugar
2742       real(kind=8) :: cou
2743       real(kind=8),dimension(3) :: ccc
2744 !el local varables
2745       integer,dimension(2,maxres/3) :: hfrag_alloc
2746       integer,dimension(4,maxres/3) :: bfrag_alloc
2747       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2748       real(kind=8),dimension(:,:), allocatable  :: c_temporary
2749       integer,dimension(:,:) , allocatable  :: itype_temporary
2750       integer,dimension(:),allocatable :: istype_temp
2751       efree_temp=0.0d0
2752       ibeg=1
2753       ishift1=0
2754       ishift=0
2755       molecule=1
2756       counter=0
2757 !      write (2,*) "UNRES_PDB",unres_pdb
2758       ires=0
2759       ires_old=0
2760       nres=0
2761       iii=0
2762       lsecondary=.false.
2763       nhfrag=0
2764       nbfrag=0
2765 !-----------------------------
2766       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2767       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2768
2769       do i=1,100000
2770         read (ipdbin,'(a80)',end=10) card
2771 !       write (iout,'(a)') card
2772         if (card(:5).eq.'HELIX') then
2773           nhfrag=nhfrag+1
2774           lsecondary=.true.
2775           read(card(22:25),*) hfrag(1,nhfrag)
2776           read(card(34:37),*) hfrag(2,nhfrag)
2777         endif
2778         if (card(:5).eq.'SHEET') then
2779           nbfrag=nbfrag+1
2780           lsecondary=.true.
2781           read(card(24:26),*) bfrag(1,nbfrag)
2782           read(card(35:37),*) bfrag(2,nbfrag)
2783 !rc----------------------------------------
2784 !rc  to be corrected !!!
2785           bfrag(3,nbfrag)=bfrag(1,nbfrag)
2786           bfrag(4,nbfrag)=bfrag(2,nbfrag)
2787 !rc----------------------------------------
2788         endif
2789         if (card(:3).eq.'END') then
2790           goto 10
2791         else if (card(:3).eq.'TER') then
2792 ! End current chain
2793           ires_old=ires+2
2794           ishift=ishift+1
2795           ishift1=ishift1+1
2796           itype(ires_old,molecule)=ntyp1_molec(molecule)
2797           itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2798           nres_molec(molecule)=nres_molec(molecule)+2
2799           ibeg=2
2800 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2801           if (unres_pdb) then
2802             do j=1,3
2803               dc(j,ires)=sccor(j,iii)
2804             enddo
2805           else
2806             call sccenter(ires,iii,sccor)
2807 !          iii=0
2808           endif
2809           iii=0
2810         endif
2811 ! Read free energy
2812         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2813 ! Fish out the ATOM cards.
2814 !        write(iout,*) 'card',card(1:20)
2815         if (index(card(1:4),'ATOM').gt.0) then
2816           read (card(12:16),*) atom
2817 !          write (iout,*) "! ",atom," !",ires
2818 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2819           read (card(23:26),*) ires
2820           read (card(18:20),'(a3)') res
2821 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2822 !     &      " ires_old",ires_old
2823 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2824 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2825           if (ires-ishift+ishift1.ne.ires_old) then
2826 ! Calculate the CM of the preceding residue.
2827 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2828             if (ibeg.eq.0) then
2829 !              write (iout,*) "Calculating sidechain center iii",iii
2830               if (unres_pdb) then
2831                 do j=1,3
2832                   dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2833                 enddo
2834               else
2835                 call sccenter(ires_old,iii,sccor)
2836               endif
2837               iii=0
2838             endif
2839 ! Start new residue.
2840             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2841               ires=ires_old
2842               cycle
2843             else if (ibeg.eq.1) then
2844               write (iout,*) "BEG ires",ires
2845               ishift=ires-1
2846               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2847                 ishift=ishift-1
2848                 itype(1,1)=ntyp1
2849                 nres_molec(molecule)=nres_molec(molecule)+1
2850               endif
2851               ires=ires-ishift+ishift1
2852               ires_old=ires
2853 !              write (iout,*) "ishift",ishift," ires",ires,&
2854 !               " ires_old",ires_old
2855               ibeg=0
2856             else if (ibeg.eq.2) then
2857 ! Start a new chain
2858               ishift=-ires_old+ires-1 !!!!!
2859               ishift1=ishift1-1    !!!!!
2860 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2861               ires=ires-ishift+ishift1
2862 !              print *,ires,ishift,ishift1
2863               ires_old=ires
2864               ibeg=0
2865             else
2866               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2867               ires=ires-ishift+ishift1
2868               ires_old=ires
2869             endif
2870 !            print *,'atom',ires,atom
2871             if (res.eq.'ACE' .or. res.eq.'NHE') then
2872               itype(ires,1)=10
2873             else
2874              if (atom.eq.'CA  '.or.atom.eq.'N   ') then
2875              molecule=1
2876               itype(ires,molecule)=rescode(ires,res,0,molecule)
2877               firstion=0
2878 !              nres_molec(molecule)=nres_molec(molecule)+1
2879             else
2880              molecule=2
2881               itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2882 !              nres_molec(molecule)=nres_molec(molecule)+1
2883              read (card(19:19),'(a1)') sugar
2884              isugar=sugarcode(sugar,ires)
2885 !            if (ibeg.eq.1) then
2886 !              istype(1)=isugar
2887 !            else
2888               istype(ires)=isugar
2889 !              print *,"ires=",ires,istype(ires)
2890 !            endif
2891
2892             endif
2893             endif
2894           else
2895             ires=ires-ishift+ishift1
2896           endif
2897 !          write (iout,*) "ires_old",ires_old," ires",ires
2898           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2899 !            ishift1=ishift1+1
2900           endif
2901 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2902           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2903              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2904             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2905 !              print *,ires,ishift,ishift1
2906 !            write (iout,*) "backbone ",atom
2907 #ifdef DEBUG
2908             write (iout,'(2i3,2x,a,3f8.3)') &
2909             ires,itype(ires,1),res,(c(j,ires),j=1,3)
2910 #endif
2911             iii=iii+1
2912               nres_molec(molecule)=nres_molec(molecule)+1
2913             do j=1,3
2914               sccor(j,iii)=c(j,ires)
2915             enddo
2916           else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2917                atom.eq."C2'" .or. atom.eq."C3'" &
2918                .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2919             read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2920 !c            write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2921 !              print *,ires,ishift,ishift1
2922             counter=counter+1
2923 !            iii=iii+1
2924 !            do j=1,3
2925 !              sccor(j,iii)=c(j,ires)
2926 !            enddo
2927             do j=1,3
2928               c(j,ires)=c(j,ires)+ccc(j)/5.0
2929             enddo
2930              print *,counter,molecule
2931              if (counter.eq.5) then
2932 !            iii=iii+1
2933               nres_molec(molecule)=nres_molec(molecule)+1
2934               firstion=0
2935 !            do j=1,3
2936 !              sccor(j,iii)=c(j,ires)
2937 !            enddo
2938              counter=0
2939            endif
2940 !            print *, "ATOM",atom(1:3)
2941           else if (atom.eq."C5'") then
2942              read (card(19:19),'(a1)') sugar
2943              isugar=sugarcode(sugar,ires)
2944             if (ibeg.eq.1) then
2945               istype(1)=isugar
2946             else
2947               istype(ires)=isugar
2948 !              print *,ires,istype(ires)
2949             endif
2950             if (unres_pdb) then
2951               read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2952             else
2953               iii=iii+1
2954               read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2955             endif
2956 !            write (*,*) card(23:27),ires,itype(ires,1)
2957           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2958                    atom.ne.'N' .and. atom.ne.'C' .and. &
2959                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2960                    atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2961                    .and. atom.ne.'P  '.and. &
2962                   atom(1:1).ne.'H' .and. &
2963                   atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2964                   ) then
2965 !            write (iout,*) "sidechain ",atom
2966 !            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2967                  if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2968 !                        write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2969             iii=iii+1
2970             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2971               endif
2972           endif
2973         else if ((ions).and.(card(1:6).eq.'HETATM')) then
2974        if (firstion.eq.0) then
2975        firstion=1
2976        if (unres_pdb) then
2977          do j=1,3
2978            dc(j,ires)=sccor(j,iii)
2979          enddo
2980        else
2981           call sccenter(ires,iii,sccor)
2982        endif
2983        endif
2984           read (card(12:16),*) atom
2985           print *,"HETATOM", atom
2986           read (card(18:20),'(a3)') res
2987           if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
2988           (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
2989           .or.(atom(1:2).eq.'K ')) &
2990           then
2991            ires=ires+1
2992            if (molecule.ne.5) molecprev=molecule
2993            molecule=5
2994            nres_molec(molecule)=nres_molec(molecule)+1
2995            print *,"HERE",nres_molec(molecule)
2996            itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2997            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2998           endif
2999         endif !atom
3000       enddo
3001    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3002       if (ires.eq.0) return
3003 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3004 ! system
3005       nres=ires
3006       if ((ires_old.ne.ires).and.(molecule.ne.5)) &
3007          nres_molec(molecule)=nres_molec(molecule)-2
3008 !      print *,'I have', nres_molec(:)
3009
3010       do k=1,4 ! ions are without dummy
3011        if (nres_molec(k).eq.0) cycle
3012       do i=2,nres-1
3013 !        write (iout,*) i,itype(i,1)
3014 !        if (itype(i,1).eq.ntyp1) then
3015 !          write (iout,*) "dummy",i,itype(i,1)
3016 !          do j=1,3
3017 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3018 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
3019 !            dc(j,i)=c(j,i)
3020 !          enddo
3021 !        endif
3022         if (itype(i,k).eq.ntyp1_molec(k)) then
3023          if (itype(i+1,k).eq.ntyp1_molec(k)) then
3024           if (itype(i+2,k).eq.0) then
3025 !           print *,"masz sieczke"
3026            do j=1,5
3027             if (itype(i+2,j).ne.0) then
3028             itype(i+1,k)=0
3029             itype(i+1,j)=ntyp1_molec(j)
3030             nres_molec(k)=nres_molec(k)-1
3031             nres_molec(j)=nres_molec(j)+1
3032             go to 3331
3033             endif
3034            enddo
3035  3331      continue
3036           endif
3037 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3038 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3039 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3040            if (unres_pdb) then
3041 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3042 !            print *,i,'tu dochodze'
3043             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3044             if (fail) then
3045               e2(1)=0.0d0
3046               e2(2)=1.0d0
3047               e2(3)=0.0d0
3048             endif !fail
3049 !            print *,i,'a tu?'
3050             do j=1,3
3051              c(j,i)=c(j,i-1)-1.9d0*e2(j)
3052             enddo
3053            else   !unres_pdb
3054            do j=1,3
3055              dcj=(c(j,i-2)-c(j,i-3))/2.0
3056             if (dcj.eq.0) dcj=1.23591524223
3057              c(j,i)=c(j,i-1)+dcj
3058              c(j,nres+i)=c(j,i)
3059            enddo
3060           endif   !unres_pdb
3061          else     !itype(i+1,1).eq.ntyp1
3062           if (unres_pdb) then
3063 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3064             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3065             if (fail) then
3066               e2(1)=0.0d0
3067               e2(2)=1.0d0
3068               e2(3)=0.0d0
3069             endif
3070             do j=1,3
3071               c(j,i)=c(j,i+1)-1.9d0*e2(j)
3072             enddo
3073           else !unres_pdb
3074            do j=1,3
3075             dcj=(c(j,i+3)-c(j,i+2))/2.0
3076             if (dcj.eq.0) dcj=1.23591524223
3077             c(j,i)=c(j,i+1)-dcj
3078             c(j,nres+i)=c(j,i)
3079            enddo
3080           endif !unres_pdb
3081          endif !itype(i+1,1).eq.ntyp1
3082         endif  !itype.eq.ntyp1
3083
3084       enddo
3085       enddo
3086 ! Calculate the CM of the last side chain.
3087       if (iii.gt.0)  then
3088       if (unres_pdb) then
3089         do j=1,3
3090           dc(j,ires)=sccor(j,iii)
3091         enddo
3092       else
3093         call sccenter(ires,iii,sccor)
3094       endif
3095       endif
3096 !      nres=ires
3097       nsup=nres
3098       nstart_sup=1
3099 !      print *,"molecule",molecule
3100       if ((itype(nres,1).ne.10)) then
3101         nres=nres+1
3102           if (molecule.eq.5) molecule=molecprev
3103         itype(nres,molecule)=ntyp1_molec(molecule)
3104         nres_molec(molecule)=nres_molec(molecule)+1
3105         if (unres_pdb) then
3106 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3107           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3108           if (fail) then
3109             e2(1)=0.0d0
3110             e2(2)=1.0d0
3111             e2(3)=0.0d0
3112           endif
3113           do j=1,3
3114             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3115           enddo
3116         else
3117         do j=1,3
3118           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3119           c(j,nres)=c(j,nres-1)+dcj
3120           c(j,2*nres)=c(j,nres)
3121         enddo
3122         endif
3123       endif
3124 !     print *,'I have',nres, nres_molec(:)
3125
3126 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3127 #ifdef WHAM_RUN
3128       if (nres.ne.nres0) then
3129         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3130                        " NRES0=",nres0
3131         stop "Error nres value in WHAM input"
3132       endif
3133 #endif
3134 !---------------------------------
3135 !el reallocate tables
3136 !      do i=1,maxres/3
3137 !       do j=1,2
3138 !         hfrag_alloc(j,i)=hfrag(j,i)
3139 !        enddo
3140 !       do j=1,4
3141 !         bfrag_alloc(j,i)=bfrag(j,i)
3142 !        enddo
3143 !      enddo
3144
3145 !      deallocate(hfrag)
3146 !      deallocate(bfrag)
3147 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
3148 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3149 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3150 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
3151
3152 !      do i=1,nhfrag
3153 !       do j=1,2
3154 !         hfrag(j,i)=hfrag_alloc(j,i)
3155 !        enddo
3156 !      enddo
3157 !      do i=1,nbfrag
3158 !       do j=1,4
3159 !         bfrag(j,i)=bfrag_alloc(j,i)
3160 !        enddo
3161 !      enddo
3162 !el end reallocate tables
3163 !---------------------------------
3164       do i=2,nres-1
3165         do j=1,3
3166           c(j,i+nres)=dc(j,i)
3167         enddo
3168       enddo
3169       do j=1,3
3170         c(j,nres+1)=c(j,1)
3171         c(j,2*nres)=c(j,nres)
3172       enddo
3173
3174       if (itype(1,1).eq.ntyp1) then
3175         nsup=nsup-1
3176         nstart_sup=2
3177         if (unres_pdb) then
3178 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3179           call refsys(2,3,4,e1,e2,e3,fail)
3180           if (fail) then
3181             e2(1)=0.0d0
3182             e2(2)=1.0d0
3183             e2(3)=0.0d0
3184           endif
3185           do j=1,3
3186             c(j,1)=c(j,2)-1.9d0*e2(j)
3187           enddo
3188         else
3189         do j=1,3
3190           dcj=(c(j,4)-c(j,3))/2.0
3191           c(j,1)=c(j,2)-dcj
3192           c(j,nres+1)=c(j,1)
3193         enddo
3194         endif
3195       endif
3196 ! First lets assign correct dummy to correct type of chain
3197 ! 1) First residue
3198       if (itype(1,1).eq.ntyp1) then
3199         if (itype(2,1).eq.0) then
3200          do j=2,5
3201            if (itype(2,j).ne.0) then
3202            itype(1,1)=0
3203            itype(1,j)=ntyp1_molec(j)
3204            nres_molec(1)=nres_molec(1)-1
3205            nres_molec(j)=nres_molec(j)+1
3206            go to 3231
3207            endif
3208          enddo
3209 3231    continue
3210         endif
3211        endif
3212        print *,'I have',nres, nres_molec(:)
3213
3214 ! Copy the coordinates to reference coordinates
3215 !      do i=1,2*nres
3216 !        do j=1,3
3217 !          cref(j,i)=c(j,i)
3218 !        enddo
3219 !      enddo
3220 ! Calculate internal coordinates.
3221       if (lprn) then
3222       write (iout,'(/a)') &
3223         "Cartesian coordinates of the reference structure"
3224       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3225        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3226       do ires=1,nres
3227         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3228           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3229           (c(j,ires+nres),j=1,3)
3230       enddo
3231       endif
3232 ! znamy już nres wiec mozna alokowac tablice
3233 ! Calculate internal coordinates.
3234       if(me.eq.king.or..not.out1file)then
3235        write (iout,'(a)') &
3236          "Backbone and SC coordinates as read from the PDB"
3237        do ires=1,nres
3238         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3239           ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3240           (c(j,nres+ires),j=1,3)
3241        enddo
3242       endif
3243 ! NOW LETS ROCK! SORTING
3244       allocate(c_temporary(3,2*nres))
3245       allocate(itype_temporary(nres,5))
3246       allocate(molnum(nres))
3247       allocate(istype_temp(nres))
3248        itype_temporary(:,:)=0
3249       seqalingbegin=1
3250       do k=1,5
3251         do i=1,nres
3252          if (itype(i,k).ne.0) then
3253           do j=1,3
3254           c_temporary(j,seqalingbegin)=c(j,i)
3255           c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3256
3257           enddo
3258           itype_temporary(seqalingbegin,k)=itype(i,k)
3259           print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3260           istype_temp(seqalingbegin)=istype(i)
3261           molnum(seqalingbegin)=k
3262           seqalingbegin=seqalingbegin+1
3263          endif
3264         enddo
3265        enddo
3266        do i=1,2*nres
3267         do j=1,3
3268         c(j,i)=c_temporary(j,i)
3269         enddo
3270        enddo
3271        do k=1,5
3272         do i=1,nres
3273          itype(i,k)=itype_temporary(i,k)
3274          istype(i)=istype_temp(i)
3275         enddo
3276        enddo
3277       if (itype(1,1).eq.ntyp1) then
3278         nsup=nsup-1
3279         nstart_sup=2
3280         if (unres_pdb) then
3281 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3282           call refsys(2,3,4,e1,e2,e3,fail)
3283           if (fail) then
3284             e2(1)=0.0d0
3285             e2(2)=1.0d0
3286             e2(3)=0.0d0
3287           endif
3288           do j=1,3
3289             c(j,1)=c(j,2)-1.9d0*e2(j)
3290           enddo
3291         else
3292         do j=1,3
3293           dcj=(c(j,4)-c(j,3))/2.0
3294           c(j,1)=c(j,2)-dcj
3295           c(j,nres+1)=c(j,1)
3296         enddo
3297         endif
3298       endif
3299
3300       if (lprn) then
3301       write (iout,'(/a)') &
3302         "Cartesian coordinates of the reference structure after sorting"
3303       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3304        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3305       do ires=1,nres
3306         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3307           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3308           (c(j,ires+nres),j=1,3)
3309       enddo
3310       endif
3311
3312 !       print *,seqalingbegin,nres
3313       if(.not.allocated(vbld)) then
3314        allocate(vbld(2*nres))
3315        do i=1,2*nres
3316          vbld(i)=0.d0
3317        enddo
3318       endif
3319       if(.not.allocated(vbld_inv)) then
3320        allocate(vbld_inv(2*nres))
3321        do i=1,2*nres
3322          vbld_inv(i)=0.d0
3323        enddo
3324       endif
3325 !!!el
3326       if(.not.allocated(theta)) then
3327         allocate(theta(nres+2))
3328         theta(:)=0.0d0
3329       endif
3330
3331       if(.not.allocated(phi)) allocate(phi(nres+2))
3332       if(.not.allocated(alph)) allocate(alph(nres+2))
3333       if(.not.allocated(omeg)) allocate(omeg(nres+2))
3334       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3335       if(.not.allocated(phiref)) allocate(phiref(nres+2))
3336       if(.not.allocated(costtab)) allocate(costtab(nres))
3337       if(.not.allocated(sinttab)) allocate(sinttab(nres))
3338       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3339       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3340       if(.not.allocated(xxref)) allocate(xxref(nres))
3341       if(.not.allocated(yyref)) allocate(yyref(nres))
3342       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3343       if(.not.allocated(dc_norm)) then
3344 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3345         allocate(dc_norm(3,0:2*nres+2))
3346         dc_norm(:,:)=0.d0
3347       endif
3348
3349       call int_from_cart(.true.,.false.)
3350       call sc_loc_geom(.false.)
3351       do i=1,nres
3352         thetaref(i)=theta(i)
3353         phiref(i)=phi(i)
3354       enddo
3355 !      do i=1,2*nres
3356 !        vbld_inv(i)=0.d0
3357 !        vbld(i)=0.d0
3358 !      enddo
3359
3360       do i=1,nres-1
3361         do j=1,3
3362           dc(j,i)=c(j,i+1)-c(j,i)
3363           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3364         enddo
3365       enddo
3366       do i=2,nres-1
3367         do j=1,3
3368           dc(j,i+nres)=c(j,i+nres)-c(j,i)
3369           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3370         enddo
3371 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3372 !        vbld_inv(i+nres)
3373       enddo
3374 !      call chainbuild
3375 ! Copy the coordinates to reference coordinates
3376 ! Splits to single chain if occurs
3377
3378 !      do i=1,2*nres
3379 !        do j=1,3
3380 !          cref(j,i,cou)=c(j,i)
3381 !        enddo
3382 !      enddo
3383 !
3384       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3385       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3386       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3387 !-----------------------------
3388       kkk=1
3389       lll=0
3390       cou=1
3391         write (iout,*) "symetr", symetr
3392       do i=1,nres
3393       lll=lll+1
3394 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3395       if (i.gt.1) then
3396       if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3397       chain_length=lll-1
3398       kkk=kkk+1
3399 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3400       lll=1
3401       endif
3402       endif
3403         do j=1,3
3404           cref(j,i,cou)=c(j,i)
3405           cref(j,i+nres,cou)=c(j,i+nres)
3406           if (i.le.nres) then
3407           chain_rep(j,lll,kkk)=c(j,i)
3408           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3409           endif
3410          enddo
3411       enddo
3412       write (iout,*) chain_length
3413       if (chain_length.eq.0) chain_length=nres
3414       do j=1,3
3415       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3416       chain_rep(j,chain_length+nres,symetr) &
3417       =chain_rep(j,chain_length+nres,1)
3418       enddo
3419 ! diagnostic
3420 !       write (iout,*) "spraw lancuchy",chain_length,symetr
3421 !       do i=1,4
3422 !         do kkk=1,chain_length
3423 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3424 !         enddo
3425 !        enddo
3426 ! enddiagnostic
3427 ! makes copy of chains
3428         write (iout,*) "symetr", symetr
3429       do j=1,3
3430       dc(j,0)=c(j,1)
3431       enddo
3432
3433       if (symetr.gt.1) then
3434        call permut(symetr)
3435        nperm=1
3436        do i=1,symetr
3437        nperm=nperm*i
3438        enddo
3439        do i=1,nperm
3440        write(iout,*) (tabperm(i,kkk),kkk=1,4)
3441        enddo
3442        do i=1,nperm
3443         cou=0
3444         do kkk=1,symetr
3445          icha=tabperm(i,kkk)
3446 !         write (iout,*) i,icha
3447          do lll=1,chain_length
3448           cou=cou+1
3449            if (cou.le.nres) then
3450            do j=1,3
3451             kupa=mod(lll,chain_length)
3452             iprzes=(kkk-1)*chain_length+lll
3453             if (kupa.eq.0) kupa=chain_length
3454 !            write (iout,*) "kupa", kupa
3455             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3456             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3457           enddo
3458           endif
3459          enddo
3460         enddo
3461        enddo
3462        endif
3463 !-koniec robienia kopii
3464 ! diag
3465       do kkk=1,nperm
3466       write (iout,*) "nowa struktura", nperm
3467       do i=1,nres
3468         write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3469       cref(2,i,kkk),&
3470       cref(3,i,kkk),cref(1,nres+i,kkk),&
3471       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3472       enddo
3473   100 format (//'              alpha-carbon coordinates       ',&
3474                 '     centroid coordinates'/ &
3475                 '       ', 6X,'X',11X,'Y',11X,'Z', &
3476                                 10X,'X',11X,'Y',11X,'Z')
3477   110 format (a,'(',i3,')',6f12.5)
3478
3479       enddo
3480 !c enddiag
3481       do j=1,nbfrag
3482         do i=1,4
3483          bfrag(i,j)=bfrag(i,j)-ishift
3484         enddo
3485       enddo
3486
3487       do j=1,nhfrag
3488         do i=1,2
3489          hfrag(i,j)=hfrag(i,j)-ishift
3490         enddo
3491       enddo
3492       ishift_pdb=ishift
3493
3494       return
3495       end subroutine readpdb
3496 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3497 !-----------------------------------------------------------------------------
3498 ! readrtns_CSA.F
3499 !-----------------------------------------------------------------------------
3500       subroutine read_control
3501 !
3502 ! Read contorl data
3503 !
3504 !      use geometry_data
3505       use comm_machsw
3506       use energy_data
3507       use control_data
3508       use compare_data
3509       use MCM_data
3510       use map_data
3511       use csa_data
3512       use MD_data
3513       use MPI_data
3514       use random, only: random_init
3515 !      implicit real*8 (a-h,o-z)
3516 !      include 'DIMENSIONS'
3517 #ifdef MP
3518       use prng, only:prng_restart
3519       include 'mpif.h'
3520       logical :: OKRandom!, prng_restart
3521       real(kind=8) :: r1
3522 #endif
3523 !      include 'COMMON.IOUNITS'
3524 !      include 'COMMON.TIME1'
3525 !      include 'COMMON.THREAD'
3526 !      include 'COMMON.SBRIDGE'
3527 !      include 'COMMON.CONTROL'
3528 !      include 'COMMON.MCM'
3529 !      include 'COMMON.MAP'
3530 !      include 'COMMON.HEADER'
3531 !      include 'COMMON.CSA'
3532 !      include 'COMMON.CHAIN'
3533 !      include 'COMMON.MUCA'
3534 !      include 'COMMON.MD'
3535 !      include 'COMMON.FFIELD'
3536 !      include 'COMMON.INTERACT'
3537 !      include 'COMMON.SETUP'
3538 !el      integer :: KDIAG,ICORFL,IXDR
3539 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3540       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3541         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
3542 !      character(len=80) :: ucase
3543       character(len=640) :: controlcard
3544
3545       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3546       integer i
3547
3548       nglob_csa=0
3549       eglob_csa=1d99
3550       nmin_csa=0
3551       read (INP,'(a)') titel
3552       call card_concat(controlcard,.true.)
3553 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3554 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3555       call reada(controlcard,'SEED',seed,0.0D0)
3556       call random_init(seed)
3557 ! Set up the time limit (caution! The time must be input in minutes!)
3558       read_cart=index(controlcard,'READ_CART').gt.0
3559       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3560       call readi(controlcard,'SYM',symetr,1)
3561       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3562       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3563       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3564       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3565       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3566       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3567       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3568       call reada(controlcard,'DRMS',drms,0.1D0)
3569       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3570        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3571        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3572        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3573        write (iout,'(a,f10.1)')'DRMS    = ',drms
3574        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3575        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3576       endif
3577       call readi(controlcard,'NZ_START',nz_start,0)
3578       call readi(controlcard,'NZ_END',nz_end,0)
3579       call readi(controlcard,'IZ_SC',iz_sc,0)
3580       timlim=60.0D0*timlim
3581       safety = 60.0d0*safety
3582       timem=timlim
3583       modecalc=0
3584       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3585 !C SHIELD keyword sets if the shielding effect of side-chains is used
3586 !C 0 denots no shielding is used all peptide are equally despite the
3587 !C solvent accesible area
3588 !C 1 the newly introduced function
3589 !C 2 reseved for further possible developement
3590       call readi(controlcard,'SHIELD',shield_mode,0)
3591 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3592         write(iout,*) "shield_mode",shield_mode
3593 !C  Varibles set size of box
3594       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3595       protein=index(controlcard,"PROTEIN").gt.0
3596       ions=index(controlcard,"IONS").gt.0
3597       nucleic=index(controlcard,"NUCLEIC").gt.0
3598       write (iout,*) "with_theta_constr ",with_theta_constr
3599       AFMlog=(index(controlcard,'AFM'))
3600       selfguide=(index(controlcard,'SELFGUIDE'))
3601       print *,'AFMlog',AFMlog,selfguide,"KUPA"
3602       call readi(controlcard,'GENCONSTR',genconstr,0)
3603       call reada(controlcard,'BOXX',boxxsize,100.0d0)
3604       call reada(controlcard,'BOXY',boxysize,100.0d0)
3605       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3606       call readi(controlcard,'TUBEMOD',tubemode,0)
3607       write (iout,*) TUBEmode,"TUBEMODE"
3608       if (TUBEmode.gt.0) then
3609        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3610        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3611        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3612        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3613        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3614        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3615        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3616        buftubebot=bordtubebot+tubebufthick
3617        buftubetop=bordtubetop-tubebufthick
3618       endif
3619
3620 ! CUTOFFF ON ELECTROSTATICS
3621       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3622       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3623       write(iout,*) "R_CUT_ELE=",r_cut_ele
3624 ! Lipidic parameters
3625       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3626       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3627       if (lipthick.gt.0.0d0) then
3628        bordliptop=(boxzsize+lipthick)/2.0
3629        bordlipbot=bordliptop-lipthick
3630       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3631       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3632       buflipbot=bordlipbot+lipbufthick
3633       bufliptop=bordliptop-lipbufthick
3634       if ((lipbufthick*2.0d0).gt.lipthick) &
3635        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3636       endif !lipthick.gt.0
3637       write(iout,*) "bordliptop=",bordliptop
3638       write(iout,*) "bordlipbot=",bordlipbot
3639       write(iout,*) "bufliptop=",bufliptop
3640       write(iout,*) "buflipbot=",buflipbot
3641       write (iout,*) "SHIELD MODE",shield_mode
3642
3643 !C-------------------------
3644       minim=(index(controlcard,'MINIMIZE').gt.0)
3645       dccart=(index(controlcard,'CART').gt.0)
3646       overlapsc=(index(controlcard,'OVERLAP').gt.0)
3647       overlapsc=.not.overlapsc
3648       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3649       searchsc=.not.searchsc
3650       sideadd=(index(controlcard,'SIDEADD').gt.0)
3651       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3652       outpdb=(index(controlcard,'PDBOUT').gt.0)
3653       outmol2=(index(controlcard,'MOL2OUT').gt.0)
3654       pdbref=(index(controlcard,'PDBREF').gt.0)
3655       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3656       indpdb=index(controlcard,'PDBSTART')
3657       extconf=(index(controlcard,'EXTCONF').gt.0)
3658       call readi(controlcard,'IPRINT',iprint,0)
3659       call readi(controlcard,'MAXGEN',maxgen,10000)
3660       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3661       call readi(controlcard,"KDIAG",kdiag,0)
3662       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3663       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3664        write (iout,*) "RESCALE_MODE",rescale_mode
3665       split_ene=index(controlcard,'SPLIT_ENE').gt.0
3666       if (index(controlcard,'REGULAR').gt.0.0D0) then
3667         call reada(controlcard,'WEIDIS',weidis,0.1D0)
3668         modecalc=1
3669         refstr=.true.
3670       endif
3671       if (index(controlcard,'CHECKGRAD').gt.0) then
3672         modecalc=5
3673         if (index(controlcard,'CART').gt.0) then
3674           icheckgrad=1
3675         elseif (index(controlcard,'CARINT').gt.0) then
3676           icheckgrad=2
3677         else
3678           icheckgrad=3
3679         endif
3680       elseif (index(controlcard,'THREAD').gt.0) then
3681         modecalc=2
3682         call readi(controlcard,'THREAD',nthread,0)
3683         if (nthread.gt.0) then
3684           call reada(controlcard,'WEIDIS',weidis,0.1D0)
3685         else
3686           if (fg_rank.eq.0) &
3687           write (iout,'(a)')'A number has to follow the THREAD keyword.'
3688           stop 'Error termination in Read_Control.'
3689         endif
3690       else if (index(controlcard,'MCMA').gt.0) then
3691         modecalc=3
3692       else if (index(controlcard,'MCEE').gt.0) then
3693         modecalc=6
3694       else if (index(controlcard,'MULTCONF').gt.0) then
3695         modecalc=4
3696       else if (index(controlcard,'MAP').gt.0) then
3697         modecalc=7
3698         call readi(controlcard,'MAP',nmap,0)
3699       else if (index(controlcard,'CSA').gt.0) then
3700         modecalc=8
3701 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3702 !rc
3703 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3704 !rc
3705 !rc        modecalc=9
3706 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3707 !fmc        modecalc=10
3708       else if (index(controlcard,'SOFTREG').gt.0) then
3709         modecalc=11
3710       else if (index(controlcard,'CHECK_BOND').gt.0) then
3711         modecalc=-1
3712       else if (index(controlcard,'TEST').gt.0) then
3713         modecalc=-2
3714       else if (index(controlcard,'MD').gt.0) then
3715         modecalc=12
3716       else if (index(controlcard,'RE ').gt.0) then
3717         modecalc=14
3718       endif
3719
3720       lmuca=index(controlcard,'MUCA').gt.0
3721       call readi(controlcard,'MUCADYN',mucadyn,0)
3722       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3723       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3724        then
3725        write (iout,*) 'MUCADYN=',mucadyn
3726        write (iout,*) 'MUCASMOOTH=',muca_smooth
3727       endif
3728
3729       iscode=index(controlcard,'ONE_LETTER')
3730       indphi=index(controlcard,'PHI')
3731       indback=index(controlcard,'BACK')
3732       iranconf=index(controlcard,'RAND_CONF')
3733       i2ndstr=index(controlcard,'USE_SEC_PRED')
3734       gradout=index(controlcard,'GRADOUT').gt.0
3735       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3736       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3737       if (me.eq.king .or. .not.out1file ) &
3738         write (iout,*) "DISTCHAINMAX",distchainmax
3739
3740       if(me.eq.king.or..not.out1file) &
3741        write (iout,'(2a)') diagmeth(kdiag),&
3742         ' routine used to diagonalize matrices.'
3743       if (shield_mode.gt.0) then
3744       pi=3.141592d0
3745 !C VSolvSphere the volume of solving sphere
3746 !C      print *,pi,"pi"
3747 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3748 !C there will be no distinction between proline peptide group and normal peptide
3749 !C group in case of shielding parameters
3750       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3751       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3752       write (iout,*) VSolvSphere,VSolvSphere_div
3753 !C long axis of side chain
3754       do i=1,ntyp
3755       long_r_sidechain(i)=vbldsc0(1,i)
3756       short_r_sidechain(i)=sigma0(i)
3757       write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3758          sigma0(i)
3759       enddo
3760       buff_shield=1.0d0
3761       endif
3762       return
3763       end subroutine read_control
3764 !-----------------------------------------------------------------------------
3765       subroutine read_REMDpar
3766 !
3767 ! Read REMD settings
3768 !
3769 !       use control
3770 !       use energy
3771 !       use geometry
3772       use REMD_data
3773       use MPI_data
3774       use control_data, only:out1file
3775 !      implicit real*8 (a-h,o-z)
3776 !      include 'DIMENSIONS'
3777 !      include 'COMMON.IOUNITS'
3778 !      include 'COMMON.TIME1'
3779 !      include 'COMMON.MD'
3780       use MD_data
3781 !el #ifndef LANG0
3782 !el      include 'COMMON.LANGEVIN'
3783 !el #else
3784 !el      include 'COMMON.LANGEVIN.lang0'
3785 !el #endif
3786 !      include 'COMMON.INTERACT'
3787 !      include 'COMMON.NAMES'
3788 !      include 'COMMON.GEO'
3789 !      include 'COMMON.REMD'
3790 !      include 'COMMON.CONTROL'
3791 !      include 'COMMON.SETUP'
3792 !      character(len=80) :: ucase
3793       character(len=320) :: controlcard
3794       character(len=3200) :: controlcard1
3795       integer :: iremd_m_total
3796 !el local variables
3797       integer :: i
3798 !     real(kind=8) :: var,ene
3799
3800       if(me.eq.king.or..not.out1file) &
3801        write (iout,*) "REMD setup"
3802
3803       call card_concat(controlcard,.true.)
3804       call readi(controlcard,"NREP",nrep,3)
3805       call readi(controlcard,"NSTEX",nstex,1000)
3806       call reada(controlcard,"RETMIN",retmin,10.0d0)
3807       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3808       mremdsync=(index(controlcard,'SYNC').gt.0)
3809       call readi(controlcard,"NSYN",i_sync_step,100)
3810       restart1file=(index(controlcard,'REST1FILE').gt.0)
3811       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3812       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3813       if(max_cache_traj_use.gt.max_cache_traj) &
3814                 max_cache_traj_use=max_cache_traj
3815       if(me.eq.king.or..not.out1file) then
3816 !d       if (traj1file) then
3817 !rc caching is in testing - NTWX is not ignored
3818 !d        write (iout,*) "NTWX value is ignored"
3819 !d        write (iout,*) "  trajectory is stored to one file by master"
3820 !d        write (iout,*) "  before exchange at NSTEX intervals"
3821 !d       endif
3822        write (iout,*) "NREP= ",nrep
3823        write (iout,*) "NSTEX= ",nstex
3824        write (iout,*) "SYNC= ",mremdsync
3825        write (iout,*) "NSYN= ",i_sync_step
3826        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3827       endif
3828       remd_tlist=.false.
3829       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3830       if (index(controlcard,'TLIST').gt.0) then
3831          remd_tlist=.true.
3832          call card_concat(controlcard1,.true.)
3833          read(controlcard1,*) (remd_t(i),i=1,nrep)
3834          if(me.eq.king.or..not.out1file) &
3835           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3836       endif
3837       remd_mlist=.false.
3838       if (index(controlcard,'MLIST').gt.0) then
3839          remd_mlist=.true.
3840          call card_concat(controlcard1,.true.)
3841          read(controlcard1,*) (remd_m(i),i=1,nrep)
3842          if(me.eq.king.or..not.out1file) then
3843           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3844           iremd_m_total=0
3845           do i=1,nrep
3846            iremd_m_total=iremd_m_total+remd_m(i)
3847           enddo
3848           write (iout,*) 'Total number of replicas ',iremd_m_total
3849          endif
3850       endif
3851       if(me.eq.king.or..not.out1file) &
3852        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3853       return
3854       end subroutine read_REMDpar
3855 !-----------------------------------------------------------------------------
3856       subroutine read_MDpar
3857 !
3858 ! Read MD settings
3859 !
3860       use control_data, only: r_cut,rlamb,out1file
3861       use energy_data
3862       use geometry_data, only: pi
3863       use MPI_data
3864 !      implicit real*8 (a-h,o-z)
3865 !      include 'DIMENSIONS'
3866 !      include 'COMMON.IOUNITS'
3867 !      include 'COMMON.TIME1'
3868 !      include 'COMMON.MD'
3869       use MD_data
3870 !el #ifndef LANG0
3871 !el      include 'COMMON.LANGEVIN'
3872 !el #else
3873 !el      include 'COMMON.LANGEVIN.lang0'
3874 !el #endif
3875 !      include 'COMMON.INTERACT'
3876 !      include 'COMMON.NAMES'
3877 !      include 'COMMON.GEO'
3878 !      include 'COMMON.SETUP'
3879 !      include 'COMMON.CONTROL'
3880 !      include 'COMMON.SPLITELE'
3881 !      character(len=80) :: ucase
3882       character(len=320) :: controlcard
3883 !el local variables
3884       integer :: i,j
3885       real(kind=8) :: eta
3886
3887       call card_concat(controlcard,.true.)
3888       call readi(controlcard,"NSTEP",n_timestep,1000000)
3889       call readi(controlcard,"NTWE",ntwe,100)
3890       call readi(controlcard,"NTWX",ntwx,1000)
3891       call reada(controlcard,"DT",d_time,1.0d-1)
3892       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3893       call reada(controlcard,"DAMAX",damax,1.0d1)
3894       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3895       call readi(controlcard,"LANG",lang,0)
3896       RESPA = index(controlcard,"RESPA") .gt. 0
3897       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3898       ntime_split0=ntime_split
3899       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3900       ntime_split0=ntime_split
3901       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3902       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3903       rest = index(controlcard,"REST").gt.0
3904       tbf = index(controlcard,"TBF").gt.0
3905       usampl = index(controlcard,"USAMPL").gt.0
3906       mdpdb = index(controlcard,"MDPDB").gt.0
3907       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3908       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3909       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3910       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3911       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3912       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3913       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3914       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3915       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3916       large = index(controlcard,"LARGE").gt.0
3917       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3918       rattle = index(controlcard,"RATTLE").gt.0
3919 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3920       nset=0
3921       if(usampl) then
3922         call read_fragments
3923       endif
3924
3925       if(me.eq.king.or..not.out1file) then
3926        write (iout,*)
3927        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3928        write (iout,*)
3929        write (iout,'(a)') "The units are:"
3930        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3931        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3932         " acceleration: angstrom/(48.9 fs)**2"
3933        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3934        write (iout,*)
3935        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3936        write (iout,'(a60,f10.5,a)') &
3937         "Initial time step of numerical integration:",d_time,&
3938         " natural units"
3939        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3940        if (RESPA) then
3941         write (iout,'(2a,i4,a)') &
3942           "A-MTS algorithm used; initial time step for fast-varying",&
3943           " short-range forces split into",ntime_split," steps."
3944         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3945          r_cut," lambda",rlamb
3946        endif
3947        write (iout,'(2a,f10.5)') &
3948         "Maximum acceleration threshold to reduce the time step",&
3949         "/increase split number:",damax
3950        write (iout,'(2a,f10.5)') &
3951         "Maximum predicted energy drift to reduce the timestep",&
3952         "/increase split number:",edriftmax
3953        write (iout,'(a60,f10.5)') &
3954        "Maximum velocity threshold to reduce velocities:",dvmax
3955        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3956        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3957        if (rattle) write (iout,'(a60)') &
3958         "Rattle algorithm used to constrain the virtual bonds"
3959       endif
3960       reset_fricmat=1000
3961       if (lang.gt.0) then
3962         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3963         call reada(controlcard,"RWAT",rwat,1.4d0)
3964         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3965         surfarea=index(controlcard,"SURFAREA").gt.0
3966         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3967         if(me.eq.king.or..not.out1file)then
3968          write (iout,'(/a,$)') "Langevin dynamics calculation"
3969          if (lang.eq.1) then
3970           write (iout,'(a/)') &
3971             " with direct integration of Langevin equations"
3972          else if (lang.eq.2) then
3973           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3974          else if (lang.eq.3) then
3975           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3976          else if (lang.eq.4) then
3977           write (iout,'(a/)') " in overdamped mode"
3978          else
3979           write (iout,'(//a,i5)') &
3980             "=========== ERROR: Unknown Langevin dynamics mode:",lang
3981           stop
3982          endif
3983          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3984          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3985          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3986          write (iout,'(a60,f10.5)') &
3987          "Scaling factor of the friction forces:",scal_fric
3988          if (surfarea) write (iout,'(2a,i10,a)') &
3989            "Friction coefficients will be scaled by solvent-accessible",&
3990            " surface area every",reset_fricmat," steps."
3991         endif
3992 ! Calculate friction coefficients and bounds of stochastic forces
3993         eta=6*pi*cPoise*etawat
3994         if(me.eq.king.or..not.out1file) &
3995          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3996           eta
3997 !        allocate(gamp
3998         do j=1,5 !types of molecules
3999         gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4000         stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4001         enddo
4002         allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4003         do j=1,5 !types of molecules
4004         do i=1,ntyp
4005           gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4006           stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4007         enddo
4008         enddo
4009
4010         if(me.eq.king.or..not.out1file)then
4011          write (iout,'(/2a/)') &
4012          "Radii of site types and friction coefficients and std's of",&
4013          " stochastic forces of fully exposed sites"
4014          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4015          do i=1,ntyp
4016           write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4017            gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4018          enddo
4019         endif
4020       else if (tbf) then
4021         if(me.eq.king.or..not.out1file)then
4022          write (iout,'(a)') "Berendsen bath calculation"
4023          write (iout,'(a60,f10.5)') "Temperature:",t_bath
4024          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4025          if (reset_moment) &
4026          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4027          count_reset_moment," steps"
4028          if (reset_vel) &
4029           write (iout,'(a,i10,a)') &
4030           "Velocities will be reset at random every",count_reset_vel,&
4031          " steps"
4032         endif
4033       else
4034         if(me.eq.king.or..not.out1file) &
4035          write (iout,'(a31)') "Microcanonical mode calculation"
4036       endif
4037       if(me.eq.king.or..not.out1file)then
4038        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4039        if (usampl) then
4040           write(iout,*) "MD running with constraints."
4041           write(iout,*) "Equilibration time ", eq_time, " mtus."
4042           write(iout,*) "Constraining ", nfrag," fragments."
4043           write(iout,*) "Length of each fragment, weight and q0:"
4044           do iset=1,nset
4045            write (iout,*) "Set of restraints #",iset
4046            do i=1,nfrag
4047               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4048                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4049            enddo
4050            write(iout,*) "constraints between ", npair, "fragments."
4051            write(iout,*) "constraint pairs, weights and q0:"
4052            do i=1,npair
4053             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4054                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4055            enddo
4056            write(iout,*) "angle constraints within ", nfrag_back,&
4057             "backbone fragments."
4058            write(iout,*) "fragment, weights:"
4059            do i=1,nfrag_back
4060             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4061                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4062                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4063            enddo
4064           enddo
4065         iset=mod(kolor,nset)+1
4066        endif
4067       endif
4068       if(me.eq.king.or..not.out1file) &
4069        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4070       return
4071       end subroutine read_MDpar
4072 !-----------------------------------------------------------------------------
4073       subroutine map_read
4074
4075       use map_data
4076 !      implicit real*8 (a-h,o-z)
4077 !      include 'DIMENSIONS'
4078 !      include 'COMMON.MAP'
4079 !      include 'COMMON.IOUNITS'
4080       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4081       character(len=80) :: mapcard      !,ucase
4082 !el local variables
4083       integer :: imap
4084 !     real(kind=8) :: var,ene
4085
4086       do imap=1,nmap
4087         read (inp,'(a)') mapcard
4088         mapcard=ucase(mapcard)
4089         if (index(mapcard,'PHI').gt.0) then
4090           kang(imap)=1
4091         else if (index(mapcard,'THE').gt.0) then
4092           kang(imap)=2
4093         else if (index(mapcard,'ALP').gt.0) then
4094           kang(imap)=3
4095         else if (index(mapcard,'OME').gt.0) then
4096           kang(imap)=4
4097         else
4098           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4099           stop 'Error - illegal variable spec in MAP card.'
4100         endif
4101         call readi (mapcard,'RES1',res1(imap),0)
4102         call readi (mapcard,'RES2',res2(imap),0)
4103         if (res1(imap).eq.0) then
4104           res1(imap)=res2(imap)
4105         else if (res2(imap).eq.0) then
4106           res2(imap)=res1(imap)
4107         endif
4108         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4109           write (iout,'(a)') &
4110           'Error - illegal definition of variable group in MAP.'
4111           stop 'Error - illegal definition of variable group in MAP.'
4112         endif
4113         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4114         call reada(mapcard,'TO',ang_to(imap),0.0D0)
4115         call readi(mapcard,'NSTEP',nstep(imap),0)
4116         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4117           write (iout,'(a)') &
4118            'Illegal boundary and/or step size specification in MAP.'
4119           stop 'Illegal boundary and/or step size specification in MAP.'
4120         endif
4121       enddo ! imap
4122       return
4123       end subroutine map_read
4124 !-----------------------------------------------------------------------------
4125       subroutine csaread
4126
4127       use control_data, only: vdisulf
4128       use csa_data
4129 !      implicit real*8 (a-h,o-z)
4130 !      include 'DIMENSIONS'
4131 !      include 'COMMON.IOUNITS'
4132 !      include 'COMMON.GEO'
4133 !      include 'COMMON.CSA'
4134 !      include 'COMMON.BANK'
4135 !      include 'COMMON.CONTROL'
4136 !      character(len=80) :: ucase
4137       character(len=620) :: mcmcard
4138 !el local variables
4139 !     integer :: ntf,ik,iw_pdb
4140 !     real(kind=8) :: var,ene
4141
4142       call card_concat(mcmcard,.true.)
4143
4144       call readi(mcmcard,'NCONF',nconf,50)
4145       call readi(mcmcard,'NADD',nadd,0)
4146       call readi(mcmcard,'JSTART',jstart,1)
4147       call readi(mcmcard,'JEND',jend,1)
4148       call readi(mcmcard,'NSTMAX',nstmax,500000)
4149       call readi(mcmcard,'N0',n0,1)
4150       call readi(mcmcard,'N1',n1,6)
4151       call readi(mcmcard,'N2',n2,4)
4152       call readi(mcmcard,'N3',n3,0)
4153       call readi(mcmcard,'N4',n4,0)
4154       call readi(mcmcard,'N5',n5,0)
4155       call readi(mcmcard,'N6',n6,10)
4156       call readi(mcmcard,'N7',n7,0)
4157       call readi(mcmcard,'N8',n8,0)
4158       call readi(mcmcard,'N9',n9,0)
4159       call readi(mcmcard,'N14',n14,0)
4160       call readi(mcmcard,'N15',n15,0)
4161       call readi(mcmcard,'N16',n16,0)
4162       call readi(mcmcard,'N17',n17,0)
4163       call readi(mcmcard,'N18',n18,0)
4164
4165       vdisulf=(index(mcmcard,'DYNSS').gt.0)
4166
4167       call readi(mcmcard,'NDIFF',ndiff,2)
4168       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4169       call readi(mcmcard,'IS1',is1,1)
4170       call readi(mcmcard,'IS2',is2,8)
4171       call readi(mcmcard,'NRAN0',nran0,4)
4172       call readi(mcmcard,'NRAN1',nran1,2)
4173       call readi(mcmcard,'IRR',irr,1)
4174       call readi(mcmcard,'NSEED',nseed,20)
4175       call readi(mcmcard,'NTOTAL',ntotal,10000)
4176       call reada(mcmcard,'CUT1',cut1,2.0d0)
4177       call reada(mcmcard,'CUT2',cut2,5.0d0)
4178       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4179       call readi(mcmcard,'ICMAX',icmax,3)
4180       call readi(mcmcard,'IRESTART',irestart,0)
4181 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
4182       ntbankm=0
4183 !!bankt
4184       call reada(mcmcard,'DELE',dele,20.0d0)
4185       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4186       call readi(mcmcard,'IREF',iref,0)
4187       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4188       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4189       call readi(mcmcard,'NCONF_IN',nconf_in,0)
4190       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4191       write (iout,*) "NCONF_IN",nconf_in
4192       return
4193       end subroutine csaread
4194 !-----------------------------------------------------------------------------
4195       subroutine mcmread
4196
4197       use mcm_data
4198       use control_data, only: MaxMoveType
4199       use MD_data
4200       use minim_data
4201 !      implicit real*8 (a-h,o-z)
4202 !      include 'DIMENSIONS'
4203 !      include 'COMMON.MCM'
4204 !      include 'COMMON.MCE'
4205 !      include 'COMMON.IOUNITS'
4206 !      character(len=80) :: ucase
4207       character(len=320) :: mcmcard
4208 !el local variables
4209       integer :: i
4210 !     real(kind=8) :: var,ene
4211
4212       call card_concat(mcmcard,.true.)
4213       call readi(mcmcard,'MAXACC',maxacc,100)
4214       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4215       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4216       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4217       call readi(mcmcard,'MAXREPM',maxrepm,200)
4218       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4219       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4220       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4221       call reada(mcmcard,'E_UP',e_up,5.0D0)
4222       call reada(mcmcard,'DELTE',delte,0.1D0)
4223       call readi(mcmcard,'NSWEEP',nsweep,5)
4224       call readi(mcmcard,'NSTEPH',nsteph,0)
4225       call readi(mcmcard,'NSTEPC',nstepc,0)
4226       call reada(mcmcard,'TMIN',tmin,298.0D0)
4227       call reada(mcmcard,'TMAX',tmax,298.0D0)
4228       call readi(mcmcard,'NWINDOW',nwindow,0)
4229       call readi(mcmcard,'PRINT_MC',print_mc,0)
4230       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4231       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4232       ent_read=(index(mcmcard,'ENT_READ').gt.0)
4233       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4234       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4235       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4236       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4237       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4238       if (nwindow.gt.0) then
4239         allocate(winstart(nwindow))     !!el (maxres)
4240         allocate(winend(nwindow))       !!el
4241         allocate(winlen(nwindow))       !!el
4242         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4243         do i=1,nwindow
4244           winlen(i)=winend(i)-winstart(i)+1
4245         enddo
4246       endif
4247       if (tmax.lt.tmin) tmax=tmin
4248       if (tmax.eq.tmin) then
4249         nstepc=0
4250         nsteph=0
4251       endif
4252       if (nstepc.gt.0 .and. nsteph.gt.0) then
4253         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4254         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4255       endif
4256       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4257 ! Probabilities of different move types
4258       sumpro_type(0)=0.0D0
4259       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4260       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4261       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4262       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
4263       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4264       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4265       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4266       do i=1,MaxMoveType
4267         print *,'i',i,' sumprotype',sumpro_type(i)
4268         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4269         print *,'i',i,' sumprotype',sumpro_type(i)
4270       enddo
4271       return
4272       end subroutine mcmread
4273 !-----------------------------------------------------------------------------
4274       subroutine read_minim
4275
4276       use minim_data
4277 !      implicit real*8 (a-h,o-z)
4278 !      include 'DIMENSIONS'
4279 !      include 'COMMON.MINIM'
4280 !      include 'COMMON.IOUNITS'
4281 !      character(len=80) :: ucase
4282       character(len=320) :: minimcard
4283 !el local variables
4284 !     integer :: ntf,ik,iw_pdb
4285 !     real(kind=8) :: var,ene
4286
4287       call card_concat(minimcard,.true.)
4288       call readi(minimcard,'MAXMIN',maxmin,2000)
4289       call readi(minimcard,'MAXFUN',maxfun,5000)
4290       call readi(minimcard,'MINMIN',minmin,maxmin)
4291       call readi(minimcard,'MINFUN',minfun,maxmin)
4292       call reada(minimcard,'TOLF',tolf,1.0D-2)
4293       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4294       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4295       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4296       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4297       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4298                'Options in energy minimization:'
4299       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4300        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4301        'MinMin:',MinMin,' MinFun:',MinFun,&
4302        ' TolF:',TolF,' RTolF:',RTolF
4303       return
4304       end subroutine read_minim
4305 !-----------------------------------------------------------------------------
4306       subroutine openunits
4307
4308       use MD_data, only: usampl
4309       use csa_data
4310       use MPI_data
4311       use control_data, only:out1file
4312       use control, only: getenv_loc
4313 !      implicit real*8 (a-h,o-z)
4314 !      include 'DIMENSIONS'
4315 #ifdef MPI
4316       include 'mpif.h'
4317       character(len=16) :: form,nodename
4318       integer :: nodelen,ierror,npos
4319 #endif
4320 !      include 'COMMON.SETUP'
4321 !      include 'COMMON.IOUNITS'
4322 !      include 'COMMON.MD'
4323 !      include 'COMMON.CONTROL'
4324       integer :: lenpre,lenpot,lentmp   !,ilen
4325 !el      external ilen
4326       character(len=3) :: out1file_text !,ucase
4327       character(len=3) :: ll
4328 !el      external ucase
4329 !el local variables
4330 !     integer :: ntf,ik,iw_pdb
4331 !     real(kind=8) :: var,ene
4332 !
4333 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4334       call getenv_loc("PREFIX",prefix)
4335       pref_orig = prefix
4336       call getenv_loc("POT",pot)
4337       call getenv_loc("DIRTMP",tmpdir)
4338       call getenv_loc("CURDIR",curdir)
4339       call getenv_loc("OUT1FILE",out1file_text)
4340 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4341       out1file_text=ucase(out1file_text)
4342       if (out1file_text(1:1).eq."Y") then
4343         out1file=.true.
4344       else
4345         out1file=fg_rank.gt.0
4346       endif
4347       lenpre=ilen(prefix)
4348       lenpot=ilen(pot)
4349       lentmp=ilen(tmpdir)
4350       if (lentmp.gt.0) then
4351           write (*,'(80(1h!))')
4352           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
4353           write (*,'(80(1h!))')
4354           write (*,*)"All output files will be on node /tmp directory."
4355 #ifdef MPI
4356         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4357         if (me.eq.king) then
4358           write (*,*) "The master node is ",nodename
4359         else if (fg_rank.eq.0) then
4360           write (*,*) "I am the CG slave node ",nodename
4361         else
4362           write (*,*) "I am the FG slave node ",nodename
4363         endif
4364 #endif
4365         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4366         lenpre = lentmp+lenpre+1
4367       endif
4368       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4369 ! Get the names and open the input files
4370 #if defined(WINIFL) || defined(WINPGI)
4371       open(1,file=pref_orig(:ilen(pref_orig))// &
4372         '.inp',status='old',readonly,shared)
4373        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4374 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4375 ! Get parameter filenames and open the parameter files.
4376       call getenv_loc('BONDPAR',bondname)
4377       open (ibond,file=bondname,status='old',readonly,shared)
4378       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4379       open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4380       call getenv_loc('THETPAR',thetname)
4381       open (ithep,file=thetname,status='old',readonly,shared)
4382       call getenv_loc('ROTPAR',rotname)
4383       open (irotam,file=rotname,status='old',readonly,shared)
4384       call getenv_loc('TORPAR',torname)
4385       open (itorp,file=torname,status='old',readonly,shared)
4386       call getenv_loc('TORDPAR',tordname)
4387       open (itordp,file=tordname,status='old',readonly,shared)
4388       call getenv_loc('FOURIER',fouriername)
4389       open (ifourier,file=fouriername,status='old',readonly,shared)
4390       call getenv_loc('ELEPAR',elename)
4391       open (ielep,file=elename,status='old',readonly,shared)
4392       call getenv_loc('SIDEPAR',sidename)
4393       open (isidep,file=sidename,status='old',readonly,shared)
4394
4395       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4396       open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4397       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4398       open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4399       call getenv_loc('TORPAR_NUCL',torname_nucl)
4400       open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4401       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4402       open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4403       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4404       open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4405
4406
4407 #elif (defined CRAY) || (defined AIX)
4408       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4409         action='read')
4410 !      print *,"Processor",myrank," opened file 1"
4411       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4412 !      print *,"Processor",myrank," opened file 9"
4413 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4414 ! Get parameter filenames and open the parameter files.
4415       call getenv_loc('BONDPAR',bondname)
4416       open (ibond,file=bondname,status='old',action='read')
4417       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4418       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4419
4420 !      print *,"Processor",myrank," opened file IBOND"
4421       call getenv_loc('THETPAR',thetname)
4422       open (ithep,file=thetname,status='old',action='read')
4423 !      print *,"Processor",myrank," opened file ITHEP"
4424       call getenv_loc('ROTPAR',rotname)
4425       open (irotam,file=rotname,status='old',action='read')
4426 !      print *,"Processor",myrank," opened file IROTAM"
4427       call getenv_loc('TORPAR',torname)
4428       open (itorp,file=torname,status='old',action='read')
4429 !      print *,"Processor",myrank," opened file ITORP"
4430       call getenv_loc('TORDPAR',tordname)
4431       open (itordp,file=tordname,status='old',action='read')
4432 !      print *,"Processor",myrank," opened file ITORDP"
4433       call getenv_loc('SCCORPAR',sccorname)
4434       open (isccor,file=sccorname,status='old',action='read')
4435 !      print *,"Processor",myrank," opened file ISCCOR"
4436       call getenv_loc('FOURIER',fouriername)
4437       open (ifourier,file=fouriername,status='old',action='read')
4438 !      print *,"Processor",myrank," opened file IFOURIER"
4439       call getenv_loc('ELEPAR',elename)
4440       open (ielep,file=elename,status='old',action='read')
4441 !      print *,"Processor",myrank," opened file IELEP"
4442       call getenv_loc('SIDEPAR',sidename)
4443       open (isidep,file=sidename,status='old',action='read')
4444
4445       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4446       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4447       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4448       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4449       call getenv_loc('TORPAR_NUCL',torname_nucl)
4450       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4451       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4452       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4453       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4454       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4455
4456       call getenv_loc('LIPTRANPAR',liptranname)
4457       open (iliptranpar,file=liptranname,status='old',action='read')
4458       call getenv_loc('TUBEPAR',tubename)
4459       open (itube,file=tubename,status='old',action='read')
4460       call getenv_loc('IONPAR',ionname)
4461       open (iion,file=ionname,status='old',action='read')
4462
4463 !      print *,"Processor",myrank," opened file ISIDEP"
4464 !      print *,"Processor",myrank," opened parameter files"
4465 #elif (defined G77)
4466       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4467       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4468 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4469 ! Get parameter filenames and open the parameter files.
4470       call getenv_loc('BONDPAR',bondname)
4471       open (ibond,file=bondname,status='old')
4472       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4473       open (ibond_nucl,file=bondname_nucl,status='old')
4474
4475       call getenv_loc('THETPAR',thetname)
4476       open (ithep,file=thetname,status='old')
4477       call getenv_loc('ROTPAR',rotname)
4478       open (irotam,file=rotname,status='old')
4479       call getenv_loc('TORPAR',torname)
4480       open (itorp,file=torname,status='old')
4481       call getenv_loc('TORDPAR',tordname)
4482       open (itordp,file=tordname,status='old')
4483       call getenv_loc('SCCORPAR',sccorname)
4484       open (isccor,file=sccorname,status='old')
4485       call getenv_loc('FOURIER',fouriername)
4486       open (ifourier,file=fouriername,status='old')
4487       call getenv_loc('ELEPAR',elename)
4488       open (ielep,file=elename,status='old')
4489       call getenv_loc('SIDEPAR',sidename)
4490       open (isidep,file=sidename,status='old')
4491
4492       open (ithep_nucl,file=thetname_nucl,status='old')
4493       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4494       open (irotam_nucl,file=rotname_nucl,status='old')
4495       call getenv_loc('TORPAR_NUCL',torname_nucl)
4496       open (itorp_nucl,file=torname_nucl,status='old')
4497       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4498       open (itordp_nucl,file=tordname_nucl,status='old')
4499       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4500       open (isidep_nucl,file=sidename_nucl,status='old')
4501
4502       call getenv_loc('LIPTRANPAR',liptranname)
4503       open (iliptranpar,file=liptranname,status='old')
4504       call getenv_loc('TUBEPAR',tubename)
4505       open (itube,file=tubename,status='old')
4506       call getenv_loc('IONPAR',ionname)
4507       open (iion,file=ionname,status='old')
4508 #else
4509       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4510         readonly)
4511        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4512 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4513 ! Get parameter filenames and open the parameter files.
4514       call getenv_loc('BONDPAR',bondname)
4515       open (ibond,file=bondname,status='old',action='read')
4516       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4517       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4518       call getenv_loc('THETPAR',thetname)
4519       open (ithep,file=thetname,status='old',action='read')
4520       call getenv_loc('ROTPAR',rotname)
4521       open (irotam,file=rotname,status='old',action='read')
4522       call getenv_loc('TORPAR',torname)
4523       open (itorp,file=torname,status='old',action='read')
4524       call getenv_loc('TORDPAR',tordname)
4525       open (itordp,file=tordname,status='old',action='read')
4526       call getenv_loc('SCCORPAR',sccorname)
4527       open (isccor,file=sccorname,status='old',action='read')
4528 #ifndef CRYST_THETA
4529       call getenv_loc('THETPARPDB',thetname_pdb)
4530       print *,"thetname_pdb ",thetname_pdb
4531       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4532       print *,ithep_pdb," opened"
4533 #endif
4534       call getenv_loc('FOURIER',fouriername)
4535       open (ifourier,file=fouriername,status='old',readonly)
4536       call getenv_loc('ELEPAR',elename)
4537       open (ielep,file=elename,status='old',readonly)
4538       call getenv_loc('SIDEPAR',sidename)
4539       open (isidep,file=sidename,status='old',readonly)
4540
4541       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4542       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4543       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4544       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4545       call getenv_loc('TORPAR_NUCL',torname_nucl)
4546       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4547       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4548       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4549       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4550       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4551       call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4552       open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4553       call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4554       open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4555       call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4556       open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4557
4558
4559       call getenv_loc('LIPTRANPAR',liptranname)
4560       open (iliptranpar,file=liptranname,status='old',action='read')
4561       call getenv_loc('TUBEPAR',tubename)
4562       open (itube,file=tubename,status='old',action='read')
4563       call getenv_loc('IONPAR',ionname)
4564       open (iion,file=ionname,status='old',action='read')
4565
4566 #ifndef CRYST_SC
4567       call getenv_loc('ROTPARPDB',rotname_pdb)
4568       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4569 #endif
4570 #endif
4571       call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4572 #if defined(WINIFL) || defined(WINPGI)
4573       open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4574 #elif (defined CRAY)  || (defined AIX)
4575       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4576 #elif (defined G77)
4577       open (iscpp_nucl,file=scpname_nucl,status='old')
4578 #else
4579       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4580 #endif
4581
4582 #ifndef OLDSCP
4583 !
4584 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4585 ! Use -DOLDSCP to use hard-coded constants instead.
4586 !
4587       call getenv_loc('SCPPAR',scpname)
4588 #if defined(WINIFL) || defined(WINPGI)
4589       open (iscpp,file=scpname,status='old',readonly,shared)
4590 #elif (defined CRAY)  || (defined AIX)
4591       open (iscpp,file=scpname,status='old',action='read')
4592 #elif (defined G77)
4593       open (iscpp,file=scpname,status='old')
4594 #else
4595       open (iscpp,file=scpname,status='old',action='read')
4596 #endif
4597 #endif
4598       call getenv_loc('PATTERN',patname)
4599 #if defined(WINIFL) || defined(WINPGI)
4600       open (icbase,file=patname,status='old',readonly,shared)
4601 #elif (defined CRAY)  || (defined AIX)
4602       open (icbase,file=patname,status='old',action='read')
4603 #elif (defined G77)
4604       open (icbase,file=patname,status='old')
4605 #else
4606       open (icbase,file=patname,status='old',action='read')
4607 #endif
4608 #ifdef MPI
4609 ! Open output file only for CG processes
4610 !      print *,"Processor",myrank," fg_rank",fg_rank
4611       if (fg_rank.eq.0) then
4612
4613       if (nodes.eq.1) then
4614         npos=3
4615       else
4616         npos = dlog10(dfloat(nodes-1))+1
4617       endif
4618       if (npos.lt.3) npos=3
4619       write (liczba,'(i1)') npos
4620       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4621         //')'
4622       write (liczba,form) me
4623       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4624         liczba(:ilen(liczba))
4625       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4626         //'.int'
4627       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4628         //'.pdb'
4629       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4630         liczba(:ilen(liczba))//'.mol2'
4631       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4632         liczba(:ilen(liczba))//'.stat'
4633       if (lentmp.gt.0) &
4634         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4635             //liczba(:ilen(liczba))//'.stat')
4636       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4637         //'.rst'
4638       if(usampl) then
4639           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4640        liczba(:ilen(liczba))//'.const'
4641       endif
4642
4643       endif
4644 #else
4645       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4646       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4647       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4648       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4649       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4650       if (lentmp.gt.0) &
4651         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4652           '.stat')
4653       rest2name=prefix(:ilen(prefix))//'.rst'
4654       if(usampl) then
4655          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4656       endif
4657 #endif
4658 #if defined(AIX) || defined(PGI)
4659       if (me.eq.king .or. .not. out1file) &
4660          open(iout,file=outname,status='unknown')
4661 #ifdef DEBUG
4662       if (fg_rank.gt.0) then
4663         write (liczba,'(i3.3)') myrank/nfgtasks
4664         write (ll,'(bz,i3.3)') fg_rank
4665         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4666          status='unknown')
4667       endif
4668 #endif
4669       if(me.eq.king) then
4670        open(igeom,file=intname,status='unknown',position='append')
4671        open(ipdb,file=pdbname,status='unknown')
4672        open(imol2,file=mol2name,status='unknown')
4673        open(istat,file=statname,status='unknown',position='append')
4674       else
4675 !1out       open(iout,file=outname,status='unknown')
4676       endif
4677 #else
4678       if (me.eq.king .or. .not.out1file) &
4679           open(iout,file=outname,status='unknown')
4680 #ifdef DEBUG
4681       if (fg_rank.gt.0) then
4682         write (liczba,'(i3.3)') myrank/nfgtasks
4683         write (ll,'(bz,i3.3)') fg_rank
4684         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4685          status='unknown')
4686       endif
4687 #endif
4688       if(me.eq.king) then
4689        open(igeom,file=intname,status='unknown',access='append')
4690        open(ipdb,file=pdbname,status='unknown')
4691        open(imol2,file=mol2name,status='unknown')
4692        open(istat,file=statname,status='unknown',access='append')
4693       else
4694 !1out       open(iout,file=outname,status='unknown')
4695       endif
4696 #endif
4697       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4698       csa_seed=prefix(:lenpre)//'.CSA.seed'
4699       csa_history=prefix(:lenpre)//'.CSA.history'
4700       csa_bank=prefix(:lenpre)//'.CSA.bank'
4701       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4702       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4703       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4704 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4705       csa_int=prefix(:lenpre)//'.int'
4706       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4707       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4708       csa_in=prefix(:lenpre)//'.CSA.in'
4709 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4710 ! Write file names
4711       if (me.eq.king)then
4712       write (iout,'(80(1h-))')
4713       write (iout,'(30x,a)') "FILE ASSIGNMENT"
4714       write (iout,'(80(1h-))')
4715       write (iout,*) "Input file                      : ",&
4716         pref_orig(:ilen(pref_orig))//'.inp'
4717       write (iout,*) "Output file                     : ",&
4718         outname(:ilen(outname))
4719       write (iout,*)
4720       write (iout,*) "Sidechain potential file        : ",&
4721         sidename(:ilen(sidename))
4722 #ifndef OLDSCP
4723       write (iout,*) "SCp potential file              : ",&
4724         scpname(:ilen(scpname))
4725 #endif
4726       write (iout,*) "Electrostatic potential file    : ",&
4727         elename(:ilen(elename))
4728       write (iout,*) "Cumulant coefficient file       : ",&
4729         fouriername(:ilen(fouriername))
4730       write (iout,*) "Torsional parameter file        : ",&
4731         torname(:ilen(torname))
4732       write (iout,*) "Double torsional parameter file : ",&
4733         tordname(:ilen(tordname))
4734       write (iout,*) "SCCOR parameter file : ",&
4735         sccorname(:ilen(sccorname))
4736       write (iout,*) "Bond & inertia constant file    : ",&
4737         bondname(:ilen(bondname))
4738       write (iout,*) "Bending parameter file          : ",&
4739         thetname(:ilen(thetname))
4740       write (iout,*) "Rotamer parameter file          : ",&
4741         rotname(:ilen(rotname))
4742 !el----
4743 #ifndef CRYST_THETA
4744       write (iout,*) "Thetpdb parameter file          : ",&
4745         thetname_pdb(:ilen(thetname_pdb))
4746 #endif
4747 !el
4748       write (iout,*) "Threading database              : ",&
4749         patname(:ilen(patname))
4750       if (lentmp.ne.0) &
4751       write (iout,*)" DIRTMP                          : ",&
4752         tmpdir(:lentmp)
4753       write (iout,'(80(1h-))')
4754       endif
4755       return
4756       end subroutine openunits
4757 !-----------------------------------------------------------------------------
4758       subroutine readrst
4759
4760       use geometry_data, only: nres,dc
4761       use MD_data
4762 !      implicit real*8 (a-h,o-z)
4763 !      include 'DIMENSIONS'
4764 !      include 'COMMON.CHAIN'
4765 !      include 'COMMON.IOUNITS'
4766 !      include 'COMMON.MD'
4767 !el local variables
4768       integer ::i,j
4769 !     real(kind=8) :: var,ene
4770
4771       open(irest2,file=rest2name,status='unknown')
4772       read(irest2,*) totT,EK,potE,totE,t_bath
4773       totTafm=totT
4774 !      do i=1,2*nres
4775 ! AL 4/17/17: Now reading d_t(0,:) too
4776       do i=0,2*nres
4777          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4778       enddo
4779 !      do i=1,2*nres
4780 ! AL 4/17/17: Now reading d_c(0,:) too
4781       do i=0,2*nres
4782          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4783       enddo
4784       if(usampl) then
4785              read (irest2,*) iset
4786       endif
4787       close(irest2)
4788       return
4789       end subroutine readrst
4790 !-----------------------------------------------------------------------------
4791       subroutine read_fragments
4792
4793       use energy_data
4794 !      use geometry
4795       use control_data, only:out1file
4796       use MD_data
4797       use MPI_data
4798 !      implicit real*8 (a-h,o-z)
4799 !      include 'DIMENSIONS'
4800 #ifdef MPI
4801       include 'mpif.h'
4802 #endif
4803 !      include 'COMMON.SETUP'
4804 !      include 'COMMON.CHAIN'
4805 !      include 'COMMON.IOUNITS'
4806 !      include 'COMMON.MD'
4807 !      include 'COMMON.CONTROL'
4808 !el local variables
4809       integer :: i
4810 !     real(kind=8) :: var,ene
4811
4812       read(inp,*) nset,nfrag,npair,nfrag_back
4813
4814 !el from module energy
4815 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
4816       if(.not.allocated(wfrag_back)) then
4817         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4818         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4819
4820         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4821         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
4822
4823         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4824         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
4825       endif
4826
4827       if(me.eq.king.or..not.out1file) &
4828        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4829         " nfrag_back",nfrag_back
4830       do iset=1,nset
4831          read(inp,*) mset(iset)
4832        do i=1,nfrag
4833          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4834            qinfrag(i,iset)
4835          if(me.eq.king.or..not.out1file) &
4836           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4837            ifrag(2,i,iset), qinfrag(i,iset)
4838        enddo
4839        do i=1,npair
4840         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4841           qinpair(i,iset)
4842         if(me.eq.king.or..not.out1file) &
4843          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4844           ipair(2,i,iset), qinpair(i,iset)
4845        enddo
4846        do i=1,nfrag_back
4847         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4848            wfrag_back(3,i,iset),&
4849            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4850         if(me.eq.king.or..not.out1file) &
4851          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4852          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4853        enddo
4854       enddo
4855       return
4856       end subroutine read_fragments
4857 !-----------------------------------------------------------------------------
4858 ! shift.F       io_csa
4859 !-----------------------------------------------------------------------------
4860       subroutine csa_read
4861
4862       use csa_data
4863 !      implicit real*8 (a-h,o-z)
4864 !      include 'DIMENSIONS'
4865 !      include 'COMMON.CSA'
4866 !      include 'COMMON.BANK'
4867 !      include 'COMMON.IOUNITS'
4868 !el local variables
4869 !     integer :: ntf,ik,iw_pdb
4870 !     real(kind=8) :: var,ene
4871
4872       open(icsa_in,file=csa_in,status="old",err=100)
4873        read(icsa_in,*) nconf
4874        read(icsa_in,*) jstart,jend
4875        read(icsa_in,*) nstmax
4876        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4877        read(icsa_in,*) nran0,nran1,irr
4878        read(icsa_in,*) nseed
4879        read(icsa_in,*) ntotal,cut1,cut2
4880        read(icsa_in,*) estop
4881        read(icsa_in,*) icmax,irestart
4882        read(icsa_in,*) ntbankm,dele,difcut
4883        read(icsa_in,*) iref,rmscut,pnccut
4884        read(icsa_in,*) ndiff
4885       close(icsa_in)
4886
4887       return
4888
4889  100  continue
4890       return
4891       end subroutine csa_read
4892 !-----------------------------------------------------------------------------
4893       subroutine initial_write
4894
4895       use csa_data
4896 !      implicit real*8 (a-h,o-z)
4897 !      include 'DIMENSIONS'
4898 !      include 'COMMON.CSA'
4899 !      include 'COMMON.BANK'
4900 !      include 'COMMON.IOUNITS'
4901 !el local variables
4902 !     integer :: ntf,ik,iw_pdb
4903 !     real(kind=8) :: var,ene
4904
4905       open(icsa_seed,file=csa_seed,status="unknown")
4906        write(icsa_seed,*) "seed"
4907       close(31)
4908 #if defined(AIX) || defined(PGI)
4909        open(icsa_history,file=csa_history,status="unknown",&
4910         position="append")
4911 #else
4912        open(icsa_history,file=csa_history,status="unknown",&
4913         access="append")
4914 #endif
4915        write(icsa_history,*) nconf
4916        write(icsa_history,*) jstart,jend
4917        write(icsa_history,*) nstmax
4918        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4919        write(icsa_history,*) nran0,nran1,irr
4920        write(icsa_history,*) nseed
4921        write(icsa_history,*) ntotal,cut1,cut2
4922        write(icsa_history,*) estop
4923        write(icsa_history,*) icmax,irestart
4924        write(icsa_history,*) ntbankm,dele,difcut
4925        write(icsa_history,*) iref,rmscut,pnccut
4926        write(icsa_history,*) ndiff
4927
4928        write(icsa_history,*)
4929       close(icsa_history)
4930
4931       open(icsa_bank1,file=csa_bank1,status="unknown")
4932        write(icsa_bank1,*) 0
4933       close(icsa_bank1)
4934
4935       return
4936       end subroutine initial_write
4937 !-----------------------------------------------------------------------------
4938       subroutine restart_write
4939
4940       use csa_data
4941 !      implicit real*8 (a-h,o-z)
4942 !      include 'DIMENSIONS'
4943 !      include 'COMMON.IOUNITS'
4944 !      include 'COMMON.CSA'
4945 !      include 'COMMON.BANK'
4946 !el local variables
4947 !     integer :: ntf,ik,iw_pdb
4948 !     real(kind=8) :: var,ene
4949
4950 #if defined(AIX) || defined(PGI)
4951        open(icsa_history,file=csa_history,position="append")
4952 #else
4953        open(icsa_history,file=csa_history,access="append")
4954 #endif
4955        write(icsa_history,*)
4956        write(icsa_history,*) "This is restart"
4957        write(icsa_history,*)
4958        write(icsa_history,*) nconf
4959        write(icsa_history,*) jstart,jend
4960        write(icsa_history,*) nstmax
4961        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4962        write(icsa_history,*) nran0,nran1,irr
4963        write(icsa_history,*) nseed
4964        write(icsa_history,*) ntotal,cut1,cut2
4965        write(icsa_history,*) estop
4966        write(icsa_history,*) icmax,irestart
4967        write(icsa_history,*) ntbankm,dele,difcut
4968        write(icsa_history,*) iref,rmscut,pnccut
4969        write(icsa_history,*) ndiff
4970        write(icsa_history,*)
4971        write(icsa_history,*) "irestart is: ", irestart
4972
4973        write(icsa_history,*)
4974       close(icsa_history)
4975
4976       return
4977       end subroutine restart_write
4978 !-----------------------------------------------------------------------------
4979 ! test.F
4980 !-----------------------------------------------------------------------------
4981       subroutine write_pdb(npdb,titelloc,ee)
4982
4983 !      implicit real*8 (a-h,o-z)
4984 !      include 'DIMENSIONS'
4985 !      include 'COMMON.IOUNITS'
4986       character(len=50) :: titelloc1
4987       character*(*) :: titelloc
4988       character(len=3) :: zahl
4989       character(len=5) :: liczba5
4990       real(kind=8) :: ee
4991       integer :: npdb   !,ilen
4992 !el      external ilen
4993 !el local variables
4994       integer :: lenpre
4995 !     real(kind=8) :: var,ene
4996
4997       titelloc1=titelloc
4998       lenpre=ilen(prefix)
4999       if (npdb.lt.1000) then
5000        call numstr(npdb,zahl)
5001        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5002       else
5003         if (npdb.lt.10000) then
5004          write(liczba5,'(i1,i4)') 0,npdb
5005         else
5006          write(liczba5,'(i5)') npdb
5007         endif
5008         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5009       endif
5010       call pdbout(ee,titelloc1,ipdb)
5011       close(ipdb)
5012       return
5013       end subroutine write_pdb
5014 !-----------------------------------------------------------------------------
5015 ! thread.F
5016 !-----------------------------------------------------------------------------
5017       subroutine write_thread_summary
5018 ! Thread the sequence through a database of known structures
5019       use control_data, only: refstr
5020 !      use geometry
5021       use energy_data, only: n_ene_comp
5022       use compare_data
5023 !      implicit real*8 (a-h,o-z)
5024 !      include 'DIMENSIONS'
5025 #ifdef MPI
5026       use MPI_data      !include 'COMMON.INFO'
5027       include 'mpif.h'
5028 #endif
5029 !      include 'COMMON.CONTROL'
5030 !      include 'COMMON.CHAIN'
5031 !      include 'COMMON.DBASE'
5032 !      include 'COMMON.INTERACT'
5033 !      include 'COMMON.VAR'
5034 !      include 'COMMON.THREAD'
5035 !      include 'COMMON.FFIELD'
5036 !      include 'COMMON.SBRIDGE'
5037 !      include 'COMMON.HEADER'
5038 !      include 'COMMON.NAMES'
5039 !      include 'COMMON.IOUNITS'
5040 !      include 'COMMON.TIME1'
5041
5042       integer,dimension(maxthread) :: ip
5043       real(kind=8),dimension(0:n_ene) :: energia
5044 !el local variables
5045       integer :: i,j,ii,jj,ipj,ik,kk,ist
5046       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5047
5048       write (iout,'(30x,a/)') &
5049        '  *********** Summary threading statistics ************'
5050       write (iout,'(a)') 'Initial energies:'
5051       write (iout,'(a4,2x,a12,14a14,3a8)') &
5052        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5053        'RMSnat','NatCONT','NNCONT','RMS'
5054 ! Energy sort patterns
5055       do i=1,nthread
5056         ip(i)=i
5057       enddo
5058       do i=1,nthread-1
5059         enet=ener(n_ene-1,ip(i))
5060         jj=i
5061         do j=i+1,nthread
5062           if (ener(n_ene-1,ip(j)).lt.enet) then
5063             jj=j
5064             enet=ener(n_ene-1,ip(j))
5065           endif
5066         enddo
5067         if (jj.ne.i) then
5068           ipj=ip(jj)
5069           ip(jj)=ip(i)
5070           ip(i)=ipj
5071         endif
5072       enddo
5073       do ik=1,nthread
5074         i=ip(ik)
5075         ii=ipatt(1,i)
5076         ist=nres_base(2,ii)+ipatt(2,i)
5077         do kk=1,n_ene_comp
5078           energia(i)=ener0(kk,i)
5079         enddo
5080         etot=ener0(n_ene_comp+1,i)
5081         rmsnat=ener0(n_ene_comp+2,i)
5082         rms=ener0(n_ene_comp+3,i)
5083         frac=ener0(n_ene_comp+4,i)
5084         frac_nn=ener0(n_ene_comp+5,i)
5085
5086         if (refstr) then
5087         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5088         i,str_nam(ii),ist+1,&
5089         (energia(print_order(kk)),kk=1,nprint_ene),&
5090         etot,rmsnat,frac,frac_nn,rms
5091         else
5092         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5093         i,str_nam(ii),ist+1,&
5094         (energia(print_order(kk)),kk=1,nprint_ene),etot
5095         endif
5096       enddo
5097       write (iout,'(//a)') 'Final energies:'
5098       write (iout,'(a4,2x,a12,17a14,3a8)') &
5099        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5100        'RMSnat','NatCONT','NNCONT','RMS'
5101       do ik=1,nthread
5102         i=ip(ik)
5103         ii=ipatt(1,i)
5104         ist=nres_base(2,ii)+ipatt(2,i)
5105         do kk=1,n_ene_comp
5106           energia(kk)=ener(kk,ik)
5107         enddo
5108         etot=ener(n_ene_comp+1,i)
5109         rmsnat=ener(n_ene_comp+2,i)
5110         rms=ener(n_ene_comp+3,i)
5111         frac=ener(n_ene_comp+4,i)
5112         frac_nn=ener(n_ene_comp+5,i)
5113         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5114         i,str_nam(ii),ist+1,&
5115         (energia(print_order(kk)),kk=1,nprint_ene),&
5116         etot,rmsnat,frac,frac_nn,rms
5117       enddo
5118       write (iout,'(/a/)') 'IEXAM array:'
5119       write (iout,'(i5)') nexcl
5120       do i=1,nexcl
5121         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5122       enddo
5123       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5124        'Max. time for threading step ',max_time_for_thread,&
5125        'Average time for threading step: ',ave_time_for_thread
5126       return
5127       end subroutine write_thread_summary
5128 !-----------------------------------------------------------------------------
5129       subroutine write_stat_thread(ithread,ipattern,ist)
5130
5131       use energy_data, only: n_ene_comp
5132       use compare_data
5133 !      implicit real*8 (a-h,o-z)
5134 !      include "DIMENSIONS"
5135 !      include "COMMON.CONTROL"
5136 !      include "COMMON.IOUNITS"
5137 !      include "COMMON.THREAD"
5138 !      include "COMMON.FFIELD"
5139 !      include "COMMON.DBASE"
5140 !      include "COMMON.NAMES"
5141       real(kind=8),dimension(0:n_ene) :: energia
5142 !el local variables
5143       integer :: ithread,ipattern,ist,i
5144       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5145
5146 #if defined(AIX) || defined(PGI)
5147       open(istat,file=statname,position='append')
5148 #else
5149       open(istat,file=statname,access='append')
5150 #endif
5151       do i=1,n_ene_comp
5152         energia(i)=ener(i,ithread)
5153       enddo
5154       etot=ener(n_ene_comp+1,ithread)
5155       rmsnat=ener(n_ene_comp+2,ithread)
5156       rms=ener(n_ene_comp+3,ithread)
5157       frac=ener(n_ene_comp+4,ithread)
5158       frac_nn=ener(n_ene_comp+5,ithread)
5159       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5160         ithread,str_nam(ipattern),ist+1,&
5161         (energia(print_order(i)),i=1,nprint_ene),&
5162         etot,rmsnat,frac,frac_nn,rms
5163       close (istat)
5164       return
5165       end subroutine write_stat_thread
5166 !-----------------------------------------------------------------------------
5167 #endif
5168 !-----------------------------------------------------------------------------
5169       end module io_config