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