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)),i,c(1,i),c(2,i),&
 482 !          c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
 483 !      enddo
 484 !  100 format (//'              alpha-carbon coordinates       ',&
 485 !                '     centroid coordinates'/ &
 486 !                '       ', 6X,'X',11X,'Y',11X,'Z',&
 487 !                                10X,'X',11X,'Y',11X,'Z')
 488 !  110 format (a,'(',i3,')',6f12.5)
 489 !      return
 490 !      end subroutine cartprint
 491 !-----------------------------------------------------------------------------
 492 ! dihed_cons.F
 493 !-----------------------------------------------------------------------------
 494       subroutine secstrp2dihc
 495
 496       use geometry_data
 497       use energy_data
 498 !      implicit real*8 (a-h,o-z)
 499 !      include 'DIMENSIONS'
 500 !      include 'COMMON.GEO'
 501 !      include 'COMMON.BOUNDS'
 502 !      include 'COMMON.CHAIN'
 503 !      include 'COMMON.TORCNSTR'
 504 !      include 'COMMON.IOUNITS'
 505 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 506 !el      COMMON/SECONDARYS/secstruc
 507       character(len=80) :: line
 508       logical :: errflag
 509 !el      external ilen
 510
 511 !el  local variables
 512       integer :: i,ii,lenpre
 513
 514       allocate(secstruc(nres))
 515
 516 !dr      call getenv_loc('SECPREDFIL',secpred)
 517       lenpre=ilen(prefix)
 518       secpred=prefix(:lenpre)//'.spred'
 519
 520 #if defined(WINIFL) || defined(WINPGI)
 521       open(isecpred,file=secpred,status='old',readonly,shared)
 522 #elif (defined CRAY) || (defined AIX)
 523       open(isecpred,file=secpred,status='old',action='read')
 524 #elif (defined G77)
 525       open(isecpred,file=secpred,status='old')
 526 #else
 527       open(isecpred,file=secpred,status='old',action='read')
 528 #endif
 529 ! read secondary structure prediction from JPRED here!
 530 !      read(isecpred,'(A80)',err=100,end=100) line
 531 !      read(line,'(f10.3)',err=110) ftors
 532        read(isecpred,'(f10.3)',err=110) ftors
 533
 534       write (iout,*) 'FTORS factor =',ftors
 535 ! initialize secstruc to any
 536        do i=1,nres
 537         secstruc(i) ='-'
 538       enddo
 539       ndih_constr=0
 540       ndih_nconstr=0
 541
 542       call read_secstr_pred(isecpred,iout,errflag)
 543       if (errflag) then
 544          write(iout,*)'There is a problem with the list of secondary-',&
 545            'structure prediction'
 546          goto 100
 547       endif
 548 ! 8/13/98 Set limits to generating the dihedral angles
 549       do i=1,nres
 550         phibound(1,i)=-pi
 551         phibound(2,i)=pi
 552       enddo
 553
 554       ii=0
 555       do i=1,nres
 556         if ( secstruc(i) .eq. 'H') then
 557 ! Helix restraints for this residue
 558            ii=ii+1
 559            idih_constr(ii)=i
 560            phi0(ii) = 45.0D0*deg2rad
 561            drange(ii)= 5.0D0*deg2rad
 562            phibound(1,i) = phi0(ii)-drange(ii)
 563            phibound(2,i) = phi0(ii)+drange(ii)
 564         else if (secstruc(i) .eq. 'E') then
 565 ! strand restraints for this residue
 566            ii=ii+1
 567            idih_constr(ii)=i
 568            phi0(ii) = 180.0D0*deg2rad
 569            drange(ii)= 5.0D0*deg2rad
 570            phibound(1,i) = phi0(ii)-drange(ii)
 571            phibound(2,i) = phi0(ii)+drange(ii)
 572         else
 573 ! no restraints for this residue
 574            ndih_nconstr=ndih_nconstr+1
 575            idih_nconstr(ndih_nconstr)=i
 576         endif
 577       enddo
 578       ndih_constr=ii
 579 !      deallocate(secstruc)
 580       return
 581 100   continue
 582       write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
 583 !      deallocate(secstruc)
 584       return
 585 110   continue
 586       write(iout,'(A20)')'Error reading FTORS'
 587 !      deallocate(secstruc)
 588       return
 589       end subroutine secstrp2dihc
 590 !-----------------------------------------------------------------------------
 591       subroutine read_secstr_pred(jin,jout,errors)
 592
 593 !      implicit real*8 (a-h,o-z)
 594 !      INCLUDE 'DIMENSIONS'
 595 !      include 'COMMON.IOUNITS'
 596 !      include 'COMMON.CHAIN'
 597 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 598 !el      COMMON/SECONDARYS/secstruc
 599 !el      EXTERNAL ILEN
 600       character(len=80) :: line,line1   !,ucase
 601       logical :: errflag,errors,blankline
 602
 603 !el  local variables
 604       integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
 605             length_of_chain
 606       errors=.false.
 607       read (jin,'(a)') line
 608       write (jout,'(2a)') '> ',line(1:78)
 609       line1=ucase(line)
 610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
 611 ! correspond to the end-groups.  ADD to the secondary structure prediction "-" for the
 612 ! end-groups in the input file "*.spred"
 613
 614       iseq=1
 615       do while (index(line1,'$END').eq.0)
 616 ! Override commented lines.
 617          ipos=1
 618          blankline=.false.
 619          do while (.not.blankline)
 620             line1=' '
 621             call mykey(line,line1,ipos,blankline,errflag)
 622             if (errflag) write (jout,'(2a)') &
 623        'Error when reading sequence in line: ',line
 624             errors=errors .or. errflag
 625             if (.not. blankline .and. .not. errflag) then
 626                ipos1=2
 627                iend=ilen(line1)
 628 !el               if (iseq.le.maxres) then
 629                   if (line1(1:1).eq.'-' ) then
 630                      secstruc(iseq)=line1(1:1)
 631                   else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
 632                             ( ucase(line1(1:1)).eq.'H' ) ) then
 633                      secstruc(iseq)=ucase(line1(1:1))
 634                   else
 635                      errors=.true.
 636                      write (jout,1010) line1(1:1), iseq
 637                      goto 80
 638                   endif
 639 !el               else
 640 !el                  errors=.true.
 641 !el                  write (jout,1000) iseq,maxres
 642 !el                  goto 80
 643 !el               endif
 644                do while (ipos1.le.iend)
 645
 646                   iseq=iseq+1
 647                   il=1
 648                   ipos1=ipos1+1
 649 !el                  if (iseq.le.maxres) then
 650                      if (line1(ipos1-1:ipos1-1).eq.'-' ) then
 651                         secstruc(iseq)=line1(ipos1-1:ipos1-1)
 652                      else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
 653                            (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
 654                         secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
 655                      else
 656                         errors=.true.
 657                         write (jout,1010) line1(ipos1-1:ipos1-1), iseq
 658                         goto 80
 659                      endif
 660 !el                  else
 661 !el                     errors=.true.
 662 !el                     write (jout,1000) iseq,maxres
 663 !el                     goto 80
 664 !el                  endif
 665                enddo
 666                iseq=iseq+1
 667             endif
 668          enddo
 669          read (jin,'(a)') line
 670          write (jout,'(2a)') '> ',line(1:78)
 671          line1=ucase(line)
 672       enddo
 673
 674 !d    write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
 675
 676 !d check whether the found length of the chain is correct.
 677       length_of_chain=iseq-1
 678       if (length_of_chain .ne. nres) then
 679 !        errors=.true.
 680         write (jout,'(a,i4,a,i4,a)') &
 681        'Error: the number of labels specified in $SEC_STRUC_PRED (' &
 682        ,length_of_chain,') does not match with the number of residues (' &
 683        ,nres,').'
 684       endif
 685    80 continue
 686
 687  1000 format('Error - the number of residues (',i4,&
 688        ') has exceeded maximum (',i4,').')
 689  1010 format ('Error - unrecognized secondary structure label',a4,&
 690        ' in position',i4)
 691       return
 692       end subroutine read_secstr_pred
 693 !#endif
 694 !-----------------------------------------------------------------------------
 695 ! parmread.F
 696 !-----------------------------------------------------------------------------
 697       subroutine parmread
 698
 699       use geometry_data
 700       use energy_data
 701       use control_data, only:maxterm !,maxtor
 702       use MD_data
 703       use MPI_data
 704 !el      use map_data
 705       use control, only: getenv_loc
 706 !
 707 ! Read the parameters of the probability distributions of the virtual-bond
 708 ! valence angles and the side chains and energy parameters.
 709 !
 710 ! Important! Energy-term weights ARE NOT read here; they are read from the
 711 ! main input file instead, because NO defaults have yet been set for these
 712 ! parameters.
 713 !
 714 !      implicit real*8 (a-h,o-z)
 715 !      include 'DIMENSIONS'
 716 #ifdef MPI
 717       include "mpif.h"
 718       integer :: IERROR
 719 #endif
 720 !      include 'COMMON.IOUNITS'
 721 !      include 'COMMON.CHAIN'
 722 !      include 'COMMON.INTERACT'
 723 !      include 'COMMON.GEO'
 724 !      include 'COMMON.LOCAL'
 725 !      include 'COMMON.TORSION'
 726 !      include 'COMMON.SCCOR'
 727 !      include 'COMMON.SCROT'
 728 !      include 'COMMON.FFIELD'
 729 !      include 'COMMON.NAMES'
 730 !      include 'COMMON.SBRIDGE'
 731 !      include 'COMMON.MD'
 732 !      include 'COMMON.SETUP'
 733       character(len=1) :: t1,t2,t3
 734       character(len=1) :: onelett(4) = (/"G","A","P","D"/)
 735       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
 736       logical :: lprint,LaTeX
 737       real(kind=8),dimension(3,3,maxlob) :: blower      !(3,3,maxlob)
 738       real(kind=8),dimension(13) :: b
 739       character(len=3) :: lancuch       !,ucase
 740 !el  local variables
 741       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
 742       integer :: maxinter,junk,kk,ii
 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(ntyp)) !(ntyp)
 780       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 781       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 782       allocate(msc(ntyp+1)) !(ntyp+1)
 783       allocate(isc(ntyp+1)) !(ntyp+1)
 784       allocate(restok(ntyp+1)) !(ntyp+1)
 785       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 786       allocate(long_r_sidechain(ntyp))
 787       allocate(short_r_sidechain(ntyp))
 788       dsc(:)=0.0d0
 789       dsc_inv(:)=0.0d0
 790
 791 #ifdef CRYST_BOND
 792       read (ibond,*) vbldp0,akp,mp,ip,pstok
 793       do i=1,ntyp
 794         nbondterm(i)=1
 795         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
 796         dsc(i) = vbldsc0(1,i)
 797         if (i.eq.10) then
 798           dsc_inv(i)=0.0D0
 799         else
 800           dsc_inv(i)=1.0D0/dsc(i)
 801         endif
 802       enddo
 803 #else
 804       read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
 805       do i=1,ntyp
 806         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
 807          j=1,nbondterm(i)),msc(i),isc(i),restok(i)
 808         dsc(i) = vbldsc0(1,i)
 809         if (i.eq.10) then
 810           dsc_inv(i)=0.0D0
 811         else
 812           dsc_inv(i)=1.0D0/dsc(i)
 813         endif
 814       enddo
 815 #endif
 816       if (lprint) then
 817         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
 818         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
 819          'inertia','Pstok'
 820         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
 821         do i=1,ntyp
 822           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
 823             vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
 824           do j=2,nbondterm(i)
 825             write (iout,'(13x,3f10.5)') &
 826               vbldsc0(j,i),aksc(j,i),abond0(j,i)
 827           enddo
 828         enddo
 829       endif
 830 !----------------------------------------------------
 831       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
 832       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))      !(-ntyp:ntyp)
 833       allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
 834       allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
 835       allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
 836       allocate(gthet(3,-ntyp:ntyp))     !(3,-ntyp:ntyp)
 837
 838       a0thet(:)=0.0D0
 839       athet(:,:,:,:)=0.0D0
 840       bthet(:,:,:,:)=0.0D0
 841       polthet(:,:)=0.0D0
 842       gthet(:,:)=0.0D0
 843       theta0(:)=0.0D0
 844       sig0(:)=0.0D0
 845       sigc0(:)=0.0D0
 846       allocate(liptranene(ntyp))
 847 !C reading lipid parameters
 848       write (iout,*) "iliptranpar",iliptranpar
 849       call flush(iout)
 850        read(iliptranpar,*) pepliptran
 851        print *,pepliptran
 852        do i=1,ntyp
 853        read(iliptranpar,*) liptranene(i)
 854         print *,liptranene(i)
 855        enddo
 856        close(iliptranpar)
 857
 858 #ifdef CRYST_THETA
 859 !
 860 ! Read the parameters of the probability distribution/energy expression
 861 ! of the virtual-bond valence angles theta
 862 !
 863       do i=1,ntyp
 864         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
 865           (bthet(j,i,1,1),j=1,2)
 866         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 867         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 868         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 869         sigc0(i)=sigc0(i)**2
 870       enddo
 871       do i=1,ntyp
 872       athet(1,i,1,-1)=athet(1,i,1,1)
 873       athet(2,i,1,-1)=athet(2,i,1,1)
 874       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 875       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 876       athet(1,i,-1,1)=-athet(1,i,1,1)
 877       athet(2,i,-1,1)=-athet(2,i,1,1)
 878       bthet(1,i,-1,1)=bthet(1,i,1,1)
 879       bthet(2,i,-1,1)=bthet(2,i,1,1)
 880       enddo
 881       do i=-ntyp,-1
 882       a0thet(i)=a0thet(-i)
 883       athet(1,i,-1,-1)=athet(1,-i,1,1)
 884       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 885       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 886       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 887       athet(1,i,-1,1)=athet(1,-i,1,1)
 888       athet(2,i,-1,1)=-athet(2,-i,1,1)
 889       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 890       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 891       athet(1,i,1,-1)=-athet(1,-i,1,1)
 892       athet(2,i,1,-1)=athet(2,-i,1,1)
 893       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 894       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 895       theta0(i)=theta0(-i)
 896       sig0(i)=sig0(-i)
 897       sigc0(i)=sigc0(-i)
 898        do j=0,3
 899         polthet(j,i)=polthet(j,-i)
 900        enddo
 901        do j=1,3
 902          gthet(j,i)=gthet(j,-i)
 903        enddo
 904       enddo
 905
 906       close (ithep)
 907       if (lprint) then
 908       if (.not.LaTeX) then
 909         write (iout,'(a)') &
 910           'Parameters of the virtual-bond valence angles:'
 911         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
 912        '    ATHETA0   ','         A1   ','        A2    ',&
 913        '        B1    ','         B2   '
 914         do i=1,ntyp
 915           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
 916               a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
 917         enddo
 918         write (iout,'(/a/9x,5a/79(1h-))') &
 919        'Parameters of the expression for sigma(theta_c):',&
 920        '     ALPH0    ','      ALPH1   ','     ALPH2    ',&
 921        '     ALPH3    ','    SIGMA0C   '
 922         do i=1,ntyp
 923           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
 924             (polthet(j,i),j=0,3),sigc0(i)
 925         enddo
 926         write (iout,'(/a/9x,5a/79(1h-))') &
 927        'Parameters of the second gaussian:',&
 928        '    THETA0    ','     SIGMA0   ','        G1    ',&
 929        '        G2    ','         G3   '
 930         do i=1,ntyp
 931           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
 932              sig0(i),(gthet(j,i),j=1,3)
 933         enddo
 934        else
 935         write (iout,'(a)') &
 936           'Parameters of the virtual-bond valence angles:'
 937         write (iout,'(/a/9x,5a/79(1h-))') &
 938        'Coefficients of expansion',&
 939        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
 940        '   b1*10^1    ','    b2*10^1   '
 941         do i=1,ntyp
 942           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
 943               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
 944               (10*bthet(j,i,1,1),j=1,2)
 945         enddo
 946         write (iout,'(/a/9x,5a/79(1h-))') &
 947        'Parameters of the expression for sigma(theta_c):',&
 948        ' alpha0       ','  alph1       ',' alph2        ',&
 949        ' alhp3        ','   sigma0c    '
 950         do i=1,ntyp
 951           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
 952             (polthet(j,i),j=0,3),sigc0(i)
 953         enddo
 954         write (iout,'(/a/9x,5a/79(1h-))') &
 955        'Parameters of the second gaussian:',&
 956        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
 957        '        G2    ','   G3*10^1    '
 958         do i=1,ntyp
 959           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
 960              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 961         enddo
 962       endif
 963       endif
 964 #else
 965 !
 966 ! Read the parameters of Utheta determined from ab initio surfaces
 967 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 968 !
 969       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
 970         ntheterm3,nsingle,ndouble
 971       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 972
 973 !----------------------------------------------------
 974       allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
 975       allocate(aa0thet(-nthetyp-1:nthetyp+1,&
 976         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 977 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 978       allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
 979         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 980 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 981 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 982       allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 983         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 984       allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 985         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 986       allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 987         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 988       allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 989         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 990 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
 991 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 992       allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 993         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 994       allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 995         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 996 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
 997 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 998
 999       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1000       do i=-ntyp1,-1
1001         ithetyp(i)=-ithetyp(-i)
1002       enddo
1003
1004       aa0thet(:,:,:,:)=0.0d0
1005       aathet(:,:,:,:,:)=0.0d0
1006       bbthet(:,:,:,:,:,:)=0.0d0
1007       ccthet(:,:,:,:,:,:)=0.0d0
1008       ddthet(:,:,:,:,:,:)=0.0d0
1009       eethet(:,:,:,:,:,:)=0.0d0
1010       ffthet(:,:,:,:,:,:,:)=0.0d0
1011       ggthet(:,:,:,:,:,:,:)=0.0d0
1012
1013 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1014       do iblock=1,2
1015 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1016 ! VAR:1=non-glicyne non-proline 2=proline
1017 ! VAR:negative values for D-aminoacid
1018       do i=0,nthetyp
1019         do j=-nthetyp,nthetyp
1020           do k=-nthetyp,nthetyp
1021             read (ithep,'(6a)',end=111,err=111) res1
1022             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1023 ! VAR: aa0thet is variable describing the average value of Foureir
1024 ! VAR: expansion series
1025 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1026 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1027 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1028             read (ithep,*,end=111,err=111) &
1029               (aathet(l,i,j,k,iblock),l=1,ntheterm)
1030             read (ithep,*,end=111,err=111) &
1031              ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1032               (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1033               (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1034               (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1035               ll=1,ntheterm2)
1036             read (ithep,*,end=111,err=111) &
1037             (((ffthet(llll,lll,ll,i,j,k,iblock),&
1038                ffthet(lll,llll,ll,i,j,k,iblock),&
1039                ggthet(llll,lll,ll,i,j,k,iblock),&
1040                ggthet(lll,llll,ll,i,j,k,iblock),&
1041                llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1042           enddo
1043         enddo
1044       enddo
1045 !
1046 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1047 ! coefficients of theta-and-gamma-dependent terms are zero.
1048 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1049 ! RECOMENTDED AFTER VERSION 3.3)
1050 !      do i=1,nthetyp
1051 !        do j=1,nthetyp
1052 !          do l=1,ntheterm
1053 !            aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1054 !            aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1055 !          enddo
1056 !          aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1057 !          aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1058 !        enddo
1059 !        do l=1,ntheterm
1060 !          aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1061 !        enddo
1062 !        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1063 !      enddo
1064 !      enddo
1065 ! AND COMMENT THE LOOPS BELOW
1066       do i=1,nthetyp
1067         do j=1,nthetyp
1068           do l=1,ntheterm
1069             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1070             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1071           enddo
1072           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1073           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1074         enddo
1075         do l=1,ntheterm
1076           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1077         enddo
1078         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1079       enddo
1080       enddo       !iblock
1081
1082 ! TILL HERE
1083 ! Substitution for D aminoacids from symmetry.
1084       do iblock=1,2
1085       do i=-nthetyp,0
1086         do j=-nthetyp,nthetyp
1087           do k=-nthetyp,nthetyp
1088            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1089            do l=1,ntheterm
1090            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1091            enddo
1092            do ll=1,ntheterm2
1093             do lll=1,nsingle
1094             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1095             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1096             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1097             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1098             enddo
1099           enddo
1100           do ll=1,ntheterm3
1101            do lll=2,ndouble
1102             do llll=1,lll-1
1103             ffthet(llll,lll,ll,i,j,k,iblock)= &
1104             ffthet(llll,lll,ll,-i,-j,-k,iblock)
1105             ffthet(lll,llll,ll,i,j,k,iblock)= &
1106             ffthet(lll,llll,ll,-i,-j,-k,iblock)
1107             ggthet(llll,lll,ll,i,j,k,iblock)= &
1108             -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1109             ggthet(lll,llll,ll,i,j,k,iblock)= &
1110             -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1111             enddo !ll
1112            enddo  !lll
1113           enddo   !llll
1114          enddo    !k
1115         enddo     !j
1116        enddo      !i
1117       enddo       !iblock
1118 !
1119 ! Control printout of the coefficients of virtual-bond-angle potentials
1120 !
1121       if (lprint) then
1122         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1123         do iblock=1,2
1124         do i=1,nthetyp+1
1125           do j=1,nthetyp+1
1126             do k=1,nthetyp+1
1127               write (iout,'(//4a)') &
1128                'Type ',onelett(i),onelett(j),onelett(k)
1129               write (iout,'(//a,10x,a)') " l","a[l]"
1130               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1131               write (iout,'(i2,1pe15.5)') &
1132                  (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1133             do l=1,ntheterm2
1134               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1135                 "b",l,"c",l,"d",l,"e",l
1136               do m=1,nsingle
1137                 write (iout,'(i2,4(1pe15.5))') m,&
1138                 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1139                 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1140               enddo
1141             enddo
1142             do l=1,ntheterm3
1143               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1144                 "f+",l,"f-",l,"g+",l,"g-",l
1145               do m=2,ndouble
1146                 do n=1,m-1
1147                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1148                     ffthet(n,m,l,i,j,k,iblock),&
1149                     ffthet(m,n,l,i,j,k,iblock),&
1150                     ggthet(n,m,l,i,j,k,iblock),&
1151                     ggthet(m,n,l,i,j,k,iblock)
1152                 enddo   !n
1153               enddo     !m
1154             enddo       !l
1155           enddo         !k
1156         enddo           !j
1157       enddo             !i
1158       enddo
1159       call flush(iout)
1160       endif
1161       write (2,*) "Start reading THETA_PDB",ithep_pdb
1162       do i=1,ntyp
1163 !      write (2,*) 'i=',i
1164         read (ithep_pdb,*,err=111,end=111) &
1165            a0thet(i),(athet(j,i,1,1),j=1,2),&
1166           (bthet(j,i,1,1),j=1,2)
1167         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1168         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1169         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1170         sigc0(i)=sigc0(i)**2
1171       enddo
1172       do i=1,ntyp
1173       athet(1,i,1,-1)=athet(1,i,1,1)
1174       athet(2,i,1,-1)=athet(2,i,1,1)
1175       bthet(1,i,1,-1)=-bthet(1,i,1,1)
1176       bthet(2,i,1,-1)=-bthet(2,i,1,1)
1177       athet(1,i,-1,1)=-athet(1,i,1,1)
1178       athet(2,i,-1,1)=-athet(2,i,1,1)
1179       bthet(1,i,-1,1)=bthet(1,i,1,1)
1180       bthet(2,i,-1,1)=bthet(2,i,1,1)
1181       enddo
1182       do i=-ntyp,-1
1183       a0thet(i)=a0thet(-i)
1184       athet(1,i,-1,-1)=athet(1,-i,1,1)
1185       athet(2,i,-1,-1)=-athet(2,-i,1,1)
1186       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1187       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1188       athet(1,i,-1,1)=athet(1,-i,1,1)
1189       athet(2,i,-1,1)=-athet(2,-i,1,1)
1190       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1191       bthet(2,i,-1,1)=bthet(2,-i,1,1)
1192       athet(1,i,1,-1)=-athet(1,-i,1,1)
1193       athet(2,i,1,-1)=athet(2,-i,1,1)
1194       bthet(1,i,1,-1)=bthet(1,-i,1,1)
1195       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1196       theta0(i)=theta0(-i)
1197       sig0(i)=sig0(-i)
1198       sigc0(i)=sigc0(-i)
1199        do j=0,3
1200         polthet(j,i)=polthet(j,-i)
1201        enddo
1202        do j=1,3
1203          gthet(j,i)=gthet(j,-i)
1204        enddo
1205       enddo
1206       write (2,*) "End reading THETA_PDB"
1207       close (ithep_pdb)
1208 #endif
1209       close(ithep)
1210
1211 !-------------------------------------------
1212       allocate(nlob(ntyp1)) !(ntyp1)
1213       allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1214       allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1215       allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1216
1217       bsc(:,:)=0.0D0
1218       nlob(:)=0
1219       nlob(:)=0
1220       dsc(:)=0.0D0
1221       censc(:,:,:)=0.0D0
1222       gaussc(:,:,:,:)=0.0D0
1223
1224 #ifdef CRYST_SC
1225 !
1226 ! Read the parameters of the probability distribution/energy expression
1227 ! of the side chains.
1228 !
1229       do i=1,ntyp
1230         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1231         if (i.eq.10) then
1232           dsc_inv(i)=0.0D0
1233         else
1234           dsc_inv(i)=1.0D0/dsc(i)
1235         endif
1236         if (i.ne.10) then
1237         do j=1,nlob(i)
1238           do k=1,3
1239             do l=1,3
1240               blower(l,k,j)=0.0D0
1241             enddo
1242           enddo
1243         enddo
1244         bsc(1,i)=0.0D0
1245         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1246           ((blower(k,l,1),l=1,k),k=1,3)
1247         censc(1,1,-i)=censc(1,1,i)
1248         censc(2,1,-i)=censc(2,1,i)
1249         censc(3,1,-i)=-censc(3,1,i)
1250         do j=2,nlob(i)
1251           read (irotam,*,end=112,err=112) bsc(j,i)
1252           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1253                                        ((blower(k,l,j),l=1,k),k=1,3)
1254         censc(1,j,-i)=censc(1,j,i)
1255         censc(2,j,-i)=censc(2,j,i)
1256         censc(3,j,-i)=-censc(3,j,i)
1257 ! BSC is amplitude of Gaussian
1258         enddo
1259         do j=1,nlob(i)
1260           do k=1,3
1261             do l=1,k
1262               akl=0.0D0
1263               do m=1,3
1264                 akl=akl+blower(k,m,j)*blower(l,m,j)
1265               enddo
1266               gaussc(k,l,j,i)=akl
1267               gaussc(l,k,j,i)=akl
1268              if (((k.eq.3).and.(l.ne.3)) &
1269               .or.((l.eq.3).and.(k.ne.3))) then
1270                 gaussc(k,l,j,-i)=-akl
1271                 gaussc(l,k,j,-i)=-akl
1272               else
1273                 gaussc(k,l,j,-i)=akl
1274                 gaussc(l,k,j,-i)=akl
1275               endif
1276             enddo
1277           enddo
1278         enddo
1279         endif
1280       enddo
1281       close (irotam)
1282       if (lprint) then
1283         write (iout,'(/a)') 'Parameters of side-chain local geometry'
1284         do i=1,ntyp
1285           nlobi=nlob(i)
1286           if (nlobi.gt.0) then
1287             if (LaTeX) then
1288               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1289                ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1290                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1291                                    'log h',(bsc(j,i),j=1,nlobi)
1292                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1293               'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1294               do k=1,3
1295                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1296                        ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1297               enddo
1298             else
1299               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1300               write (iout,'(a,f10.4,4(16x,f10.4))') &
1301                                    'Center  ',(bsc(j,i),j=1,nlobi)
1302               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1303                  j=1,nlobi)
1304               write (iout,'(a)')
1305             endif
1306           endif
1307         enddo
1308       endif
1309 #else
1310 !
1311 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1312 ! added by Urszula Kozlowska 07/11/2007
1313 !
1314 !el Maximum number of SC local term fitting function coefficiants
1315 !el      integer,parameter :: maxsccoef=65
1316
1317       allocate(sc_parmin(65,ntyp))      !(maxsccoef,ntyp)
1318
1319       do i=1,ntyp
1320         read (irotam,*,end=112,err=112)
1321        if (i.eq.10) then
1322          read (irotam,*,end=112,err=112)
1323        else
1324          do j=1,65
1325            read(irotam,*,end=112,err=112) sc_parmin(j,i)
1326          enddo
1327        endif
1328       enddo
1329 !
1330 ! Read the parameters of the probability distribution/energy expression
1331 ! of the side chains.
1332 !
1333       write (2,*) "Start reading ROTAM_PDB"
1334       do i=1,ntyp
1335         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1336         if (i.eq.10) then
1337           dsc_inv(i)=0.0D0
1338         else
1339           dsc_inv(i)=1.0D0/dsc(i)
1340         endif
1341         if (i.ne.10) then
1342         do j=1,nlob(i)
1343           do k=1,3
1344             do l=1,3
1345               blower(l,k,j)=0.0D0
1346             enddo
1347           enddo
1348         enddo
1349         bsc(1,i)=0.0D0
1350         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1351           ((blower(k,l,1),l=1,k),k=1,3)
1352         do j=2,nlob(i)
1353           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1354           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1355                                        ((blower(k,l,j),l=1,k),k=1,3)
1356         enddo
1357         do j=1,nlob(i)
1358           do k=1,3
1359             do l=1,k
1360               akl=0.0D0
1361               do m=1,3
1362                 akl=akl+blower(k,m,j)*blower(l,m,j)
1363               enddo
1364               gaussc(k,l,j,i)=akl
1365               gaussc(l,k,j,i)=akl
1366             enddo
1367           enddo
1368         enddo
1369         endif
1370       enddo
1371       close (irotam_pdb)
1372       write (2,*) "End reading ROTAM_PDB"
1373 #endif
1374       close(irotam)
1375
1376 #ifdef CRYST_TOR
1377 !
1378 ! Read torsional parameters in old format
1379 !
1380       allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1381
1382       read (itorp,*,end=113,err=113) ntortyp,nterm_old
1383       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1384       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1385
1386 !el from energy module--------
1387       allocate(v1(nterm_old,ntortyp,ntortyp))
1388       allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1389 !el---------------------------
1390       do i=1,ntortyp
1391         do j=1,ntortyp
1392           read (itorp,'(a)')
1393           do k=1,nterm_old
1394             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1395           enddo
1396         enddo
1397       enddo
1398       close (itorp)
1399       if (lprint) then
1400         write (iout,'(/a/)') 'Torsional constants:'
1401         do i=1,ntortyp
1402           do j=1,ntortyp
1403             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1404             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1405           enddo
1406         enddo
1407       endif
1408 #else
1409 !
1410 ! Read torsional parameters
1411 !
1412       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1413       read (itorp,*,end=113,err=113) ntortyp
1414 !el from energy module---------
1415       allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1416       allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1417
1418       allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1419       allocate(vlor2(maxlor,ntortyp,ntortyp))
1420       allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1421       allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1422
1423       allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1424       allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1425 !el---------------------------
1426       nterm(:,:,:)=0
1427       nlor(:,:,:)=0
1428 !el---------------------------
1429
1430       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1431       do i=-ntyp,-1
1432        itortyp(i)=-itortyp(-i)
1433       enddo
1434       itortyp(ntyp1)=ntortyp
1435       itortyp(-ntyp1)=-ntortyp
1436       do iblock=1,2
1437       write (iout,*) 'ntortyp',ntortyp
1438       do i=0,ntortyp-1
1439         do j=-ntortyp+1,ntortyp-1
1440           read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1441                 nlor(i,j,iblock)
1442           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1443           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1444           v0ij=0.0d0
1445           si=-1.0d0
1446           do k=1,nterm(i,j,iblock)
1447             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1448             v2(k,i,j,iblock)
1449             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1450             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1451             v0ij=v0ij+si*v1(k,i,j,iblock)
1452             si=-si
1453 !         write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1454 !         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1455 !      v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1456           enddo
1457           do k=1,nlor(i,j,iblock)
1458             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1459               vlor2(k,i,j),vlor3(k,i,j)
1460             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1461           enddo
1462           v0(i,j,iblock)=v0ij
1463           v0(-i,-j,iblock)=v0ij
1464         enddo
1465       enddo
1466       enddo
1467       close (itorp)
1468       if (lprint) then
1469         write (iout,'(/a/)') 'Torsional constants:'
1470         do iblock=1,2
1471         do i=-ntortyp,ntortyp
1472           do j=-ntortyp,ntortyp
1473             write (iout,*) 'ityp',i,' jtyp',j
1474             write (iout,*) 'Fourier constants'
1475             do k=1,nterm(i,j,iblock)
1476               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1477               v2(k,i,j,iblock)
1478             enddo
1479             write (iout,*) 'Lorenz constants'
1480             do k=1,nlor(i,j,iblock)
1481               write (iout,'(3(1pe15.5))') &
1482                vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1483             enddo
1484           enddo
1485         enddo
1486         enddo
1487       endif
1488 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1489 !
1490 ! 6/23/01 Read parameters for double torsionals
1491 !
1492 !el from energy module------------
1493       allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1494       allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1495 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1496       allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1497       allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1498         !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1499       allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1500       allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1501         !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1502 !---------------------------------
1503
1504       do iblock=1,2
1505       do i=0,ntortyp-1
1506         do j=-ntortyp+1,ntortyp-1
1507           do k=-ntortyp+1,ntortyp-1
1508             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1509 !              write (iout,*) "OK onelett",
1510 !     &         i,j,k,t1,t2,t3
1511
1512             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1513               .or. t3.ne.toronelet(k)) then
1514              write (iout,*) "Error in double torsional parameter file",&
1515                i,j,k,t1,t2,t3
1516 #ifdef MPI
1517               call MPI_Finalize(Ierror)
1518 #endif
1519                stop "Error in double torsional parameter file"
1520             endif
1521            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1522                ntermd_2(i,j,k,iblock)
1523             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1524             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1525             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1526                ntermd_1(i,j,k,iblock))
1527             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1528                ntermd_1(i,j,k,iblock))
1529             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1530                ntermd_1(i,j,k,iblock))
1531             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1532                ntermd_1(i,j,k,iblock))
1533 ! Martix of D parameters for one dimesional foureir series
1534             do l=1,ntermd_1(i,j,k,iblock)
1535              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1536              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1537              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1538              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1539 !            write(iout,*) "whcodze" ,
1540 !     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1541             enddo
1542             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1543                v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1544                v2s(m,l,i,j,k,iblock),&
1545                m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1546 ! Martix of D parameters for two dimesional fourier series
1547             do l=1,ntermd_2(i,j,k,iblock)
1548              do m=1,l-1
1549              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1550              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1551              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1552              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1553              enddo!m
1554             enddo!l
1555           enddo!k
1556         enddo!j
1557       enddo!i
1558       enddo!iblock
1559       if (lprint) then
1560       write (iout,*)
1561       write (iout,*) 'Constants for double torsionals'
1562       do iblock=1,2
1563       do i=0,ntortyp-1
1564         do j=-ntortyp+1,ntortyp-1
1565           do k=-ntortyp+1,ntortyp-1
1566             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1567               ' nsingle',ntermd_1(i,j,k,iblock),&
1568               ' ndouble',ntermd_2(i,j,k,iblock)
1569             write (iout,*)
1570             write (iout,*) 'Single angles:'
1571             do l=1,ntermd_1(i,j,k,iblock)
1572               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1573                  v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1574                  v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1575                  v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1576             enddo
1577             write (iout,*)
1578             write (iout,*) 'Pairs of angles:'
1579             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1580             do l=1,ntermd_2(i,j,k,iblock)
1581               write (iout,'(i5,20f10.5)') &
1582                l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1583             enddo
1584             write (iout,*)
1585             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1586             do l=1,ntermd_2(i,j,k,iblock)
1587               write (iout,'(i5,20f10.5)') &
1588                l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1589                (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1590             enddo
1591             write (iout,*)
1592           enddo
1593         enddo
1594       enddo
1595       enddo
1596       endif
1597 #endif
1598 ! Read of Side-chain backbone correlation parameters
1599 ! Modified 11 May 2012 by Adasko
1600 !CC
1601 !
1602       read (isccor,*,end=119,err=119) nsccortyp
1603
1604 !el from module energy-------------
1605       allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1606       allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1607       allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1608       allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp))   !(maxterm_sccor,20,20)
1609 !-----------------------------------
1610 #ifdef SCCORPDB
1611 !el from module energy-------------
1612       allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1613
1614       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1615       do i=-ntyp,-1
1616         isccortyp(i)=-isccortyp(-i)
1617       enddo
1618       iscprol=isccortyp(20)
1619 !      write (iout,*) 'ntortyp',ntortyp
1620       maxinter=3
1621 !c maxinter is maximum interaction sites
1622 !el from module energy---------
1623       allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1624       allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1625                -nsccortyp:nsccortyp))
1626       allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1627                -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1628       allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1629                -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1630 !-----------------------------------
1631       do l=1,maxinter
1632       do i=1,nsccortyp
1633         do j=1,nsccortyp
1634           read (isccor,*,end=119,err=119) &
1635       nterm_sccor(i,j),nlor_sccor(i,j)
1636           v0ijsccor=0.0d0
1637           v0ijsccor1=0.0d0
1638           v0ijsccor2=0.0d0
1639           v0ijsccor3=0.0d0
1640           si=-1.0d0
1641           nterm_sccor(-i,j)=nterm_sccor(i,j)
1642           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1643           nterm_sccor(i,-j)=nterm_sccor(i,j)
1644           do k=1,nterm_sccor(i,j)
1645             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1646            v2sccor(k,l,i,j)
1647             if (j.eq.iscprol) then
1648              if (i.eq.isccortyp(10)) then
1649              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1650              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1651              else
1652              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1653                               +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1654              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1655                               +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1656              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1657              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1658              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1659              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1660              endif
1661             else
1662              if (i.eq.isccortyp(10)) then
1663              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1664              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1665              else
1666                if (j.eq.isccortyp(10)) then
1667              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1668              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1669                else
1670              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1671              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1672              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1673              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1674              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1675              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1676                 endif
1677                endif
1678             endif
1679             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1680             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1681             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1682             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1683             si=-si
1684           enddo
1685           do k=1,nlor_sccor(i,j)
1686             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1687               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1688             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1689       (1+vlor3sccor(k,i,j)**2)
1690           enddo
1691           v0sccor(l,i,j)=v0ijsccor
1692           v0sccor(l,-i,j)=v0ijsccor1
1693           v0sccor(l,i,-j)=v0ijsccor2
1694           v0sccor(l,-i,-j)=v0ijsccor3
1695         enddo
1696       enddo
1697       enddo
1698       close (isccor)
1699 #else
1700 !el from module energy-------------
1701       allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1702
1703       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1704 !      write (iout,*) 'ntortyp',ntortyp
1705       maxinter=3
1706 !c maxinter is maximum interaction sites
1707 !el from module energy---------
1708       allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1709       allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1710       allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1711       allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1712 !-----------------------------------
1713       do l=1,maxinter
1714       do i=1,nsccortyp
1715         do j=1,nsccortyp
1716           read (isccor,*,end=119,err=119) &
1717        nterm_sccor(i,j),nlor_sccor(i,j)
1718           v0ijsccor=0.0d0
1719           si=-1.0d0
1720
1721           do k=1,nterm_sccor(i,j)
1722             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1723            v2sccor(k,l,i,j)
1724             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1725             si=-si
1726           enddo
1727           do k=1,nlor_sccor(i,j)
1728             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1729               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1730             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1731       (1+vlor3sccor(k,i,j)**2)
1732           enddo
1733           v0sccor(l,i,j)=v0ijsccor !el ,iblock
1734         enddo
1735       enddo
1736       enddo
1737       close (isccor)
1738
1739 #endif
1740       if (lprint) then
1741         write (iout,'(/a/)') 'Torsional constants:'
1742         do i=1,nsccortyp
1743           do j=1,nsccortyp
1744             write (iout,*) 'ityp',i,' jtyp',j
1745             write (iout,*) 'Fourier constants'
1746             do k=1,nterm_sccor(i,j)
1747       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1748             enddo
1749             write (iout,*) 'Lorenz constants'
1750             do k=1,nlor_sccor(i,j)
1751               write (iout,'(3(1pe15.5))') &
1752                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1753             enddo
1754           enddo
1755         enddo
1756       endif
1757
1758 !
1759 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1760 !         interaction energy of the Gly, Ala, and Pro prototypes.
1761 !
1762
1763       if (lprint) then
1764         write (iout,*)
1765         write (iout,*) "Coefficients of the cumulants"
1766       endif
1767       read (ifourier,*) nloctyp
1768 !write(iout,*) "nloctyp",nloctyp
1769 !el from module energy-------
1770       allocate(b1(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1771       allocate(b2(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1772       allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1773       allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1774       allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1775       allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1776       allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1777       allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1778 ! el
1779         b1(1,:)=0.0d0
1780         b1(2,:)=0.0d0
1781 !--------------------------------
1782
1783       do i=0,nloctyp-1
1784         read (ifourier,*,end=115,err=115)
1785         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1786         if (lprint) then
1787           write (iout,*) 'Type',i
1788           write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1789         endif
1790         B1(1,i)  = b(3)
1791         B1(2,i)  = b(5)
1792         B1(1,-i) = b(3)
1793         B1(2,-i) = -b(5)
1794 !        b1(1,i)=0.0d0
1795 !        b1(2,i)=0.0d0
1796         B1tilde(1,i) = b(3)
1797         B1tilde(2,i) =-b(5)
1798         B1tilde(1,-i) =-b(3)
1799         B1tilde(2,-i) =b(5)
1800 !        b1tilde(1,i)=0.0d0
1801 !        b1tilde(2,i)=0.0d0
1802         B2(1,i)  = b(2)
1803         B2(2,i)  = b(4)
1804         B2(1,-i)  =b(2)
1805         B2(2,-i)  =-b(4)
1806
1807 !        b2(1,i)=0.0d0
1808 !        b2(2,i)=0.0d0
1809         CC(1,1,i)= b(7)
1810         CC(2,2,i)=-b(7)
1811         CC(2,1,i)= b(9)
1812         CC(1,2,i)= b(9)
1813         CC(1,1,-i)= b(7)
1814         CC(2,2,-i)=-b(7)
1815         CC(2,1,-i)=-b(9)
1816         CC(1,2,-i)=-b(9)
1817 !        CC(1,1,i)=0.0d0
1818 !        CC(2,2,i)=0.0d0
1819 !        CC(2,1,i)=0.0d0
1820 !        CC(1,2,i)=0.0d0
1821         Ctilde(1,1,i)=b(7)
1822         Ctilde(1,2,i)=b(9)
1823         Ctilde(2,1,i)=-b(9)
1824         Ctilde(2,2,i)=b(7)
1825         Ctilde(1,1,-i)=b(7)
1826         Ctilde(1,2,-i)=-b(9)
1827         Ctilde(2,1,-i)=b(9)
1828         Ctilde(2,2,-i)=b(7)
1829
1830 !        Ctilde(1,1,i)=0.0d0
1831 !        Ctilde(1,2,i)=0.0d0
1832 !        Ctilde(2,1,i)=0.0d0
1833 !        Ctilde(2,2,i)=0.0d0
1834         DD(1,1,i)= b(6)
1835         DD(2,2,i)=-b(6)
1836         DD(2,1,i)= b(8)
1837         DD(1,2,i)= b(8)
1838         DD(1,1,-i)= b(6)
1839         DD(2,2,-i)=-b(6)
1840         DD(2,1,-i)=-b(8)
1841         DD(1,2,-i)=-b(8)
1842 !        DD(1,1,i)=0.0d0
1843 !        DD(2,2,i)=0.0d0
1844 !        DD(2,1,i)=0.0d0
1845 !        DD(1,2,i)=0.0d0
1846         Dtilde(1,1,i)=b(6)
1847         Dtilde(1,2,i)=b(8)
1848         Dtilde(2,1,i)=-b(8)
1849         Dtilde(2,2,i)=b(6)
1850         Dtilde(1,1,-i)=b(6)
1851         Dtilde(1,2,-i)=-b(8)
1852         Dtilde(2,1,-i)=b(8)
1853         Dtilde(2,2,-i)=b(6)
1854
1855 !        Dtilde(1,1,i)=0.0d0
1856 !        Dtilde(1,2,i)=0.0d0
1857 !        Dtilde(2,1,i)=0.0d0
1858 !        Dtilde(2,2,i)=0.0d0
1859         EE(1,1,i)= b(10)+b(11)
1860         EE(2,2,i)=-b(10)+b(11)
1861         EE(2,1,i)= b(12)-b(13)
1862         EE(1,2,i)= b(12)+b(13)
1863         EE(1,1,-i)= b(10)+b(11)
1864         EE(2,2,-i)=-b(10)+b(11)
1865         EE(2,1,-i)=-b(12)+b(13)
1866         EE(1,2,-i)=-b(12)-b(13)
1867
1868 !        ee(1,1,i)=1.0d0
1869 !        ee(2,2,i)=1.0d0
1870 !        ee(2,1,i)=0.0d0
1871 !        ee(1,2,i)=0.0d0
1872 !        ee(2,1,i)=ee(1,2,i)
1873       enddo
1874       if (lprint) then
1875       do i=1,nloctyp
1876         write (iout,*) 'Type',i
1877         write (iout,*) 'B1'
1878         write(iout,*) B1(1,i),B1(2,i)
1879         write (iout,*) 'B2'
1880         write(iout,*) B2(1,i),B2(2,i)
1881         write (iout,*) 'CC'
1882         do j=1,2
1883           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1884         enddo
1885         write(iout,*) 'DD'
1886         do j=1,2
1887           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1888         enddo
1889         write(iout,*) 'EE'
1890         do j=1,2
1891           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1892         enddo
1893       enddo
1894       endif
1895 !
1896 ! Read electrostatic-interaction parameters
1897 !
1898
1899       if (lprint) then
1900         write (iout,*)
1901         write (iout,'(/a)') 'Electrostatic interaction constants:'
1902         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1903                   'IT','JT','APP','BPP','AEL6','AEL3'
1904       endif
1905       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1906       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1907       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1908       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1909       close (ielep)
1910       do i=1,2
1911         do j=1,2
1912         rri=rpp(i,j)**6
1913         app (i,j)=epp(i,j)*rri*rri
1914         bpp (i,j)=-2.0D0*epp(i,j)*rri
1915         ael6(i,j)=elpp6(i,j)*4.2D0**6
1916         ael3(i,j)=elpp3(i,j)*4.2D0**3
1917 !        lprint=.true.
1918         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1919                           ael6(i,j),ael3(i,j)
1920 !        lprint=.false.
1921         enddo
1922       enddo
1923 !
1924 ! Read side-chain interaction parameters.
1925 !
1926 !el from module energy - COMMON.INTERACT-------
1927       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1928       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1929       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1930       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1931       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1932       allocate(epslip(ntyp,ntyp))
1933       augm(:,:)=0.0D0
1934       chip(:)=0.0D0
1935       alp(:)=0.0D0
1936       sigma0(:)=0.0D0
1937       sigii(:)=0.0D0
1938       rr0(:)=0.0D0
1939
1940 !--------------------------------
1941
1942       read (isidep,*,end=117,err=117) ipot,expon
1943       if (ipot.lt.1 .or. ipot.gt.5) then
1944         write (iout,'(2a)') 'Error while reading SC interaction',&
1945                      'potential file - unknown potential type.'
1946 #ifdef MPI
1947         call MPI_Finalize(Ierror)
1948 #endif
1949         stop
1950       endif
1951       expon2=expon/2
1952       if(me.eq.king) &
1953        write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1954        ', exponents are ',expon,2*expon
1955 !      goto (10,20,30,30,40) ipot
1956       select case(ipot)
1957 !----------------------- LJ potential ---------------------------------
1958        case (1)
1959 !   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1960          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1961            (sigma0(i),i=1,ntyp)
1962         if (lprint) then
1963           write (iout,'(/a/)') 'Parameters of the LJ potential:'
1964           write (iout,'(a/)') 'The epsilon array:'
1965           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1966           write (iout,'(/a)') 'One-body parameters:'
1967           write (iout,'(a,4x,a)') 'residue','sigma'
1968           write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1969         endif
1970 !      goto 50
1971 !----------------------- LJK potential --------------------------------
1972        case(2)
1973 !   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1974          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1975           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1976         if (lprint) then
1977           write (iout,'(/a/)') 'Parameters of the LJK potential:'
1978           write (iout,'(a/)') 'The epsilon array:'
1979           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1980           write (iout,'(/a)') 'One-body parameters:'
1981           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1982           write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1983                 i=1,ntyp)
1984         endif
1985 !      goto 50
1986 !---------------------- GB or BP potential -----------------------------
1987        case(3:4)
1988 !   30 do i=1,ntyp
1989         do i=1,ntyp
1990          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1991         enddo
1992         read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1993         read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1994         read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1995         read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1996         do i=1,ntyp
1997          read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1998         enddo
1999
2000 ! For the GB potential convert sigma'**2 into chi'
2001         if (ipot.eq.4) then
2002           do i=1,ntyp
2003             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2004           enddo
2005         endif
2006         if (lprint) then
2007           write (iout,'(/a/)') 'Parameters of the BP potential:'
2008           write (iout,'(a/)') 'The epsilon array:'
2009           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2010           write (iout,'(/a)') 'One-body parameters:'
2011           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
2012                '    chip  ','    alph  '
2013           write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
2014                              chip(i),alp(i),i=1,ntyp)
2015         endif
2016 !      goto 50
2017 !--------------------- GBV potential -----------------------------------
2018        case(5)
2019 !   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2020         read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2021           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2022           (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2023         if (lprint) then
2024           write (iout,'(/a/)') 'Parameters of the GBV potential:'
2025           write (iout,'(a/)') 'The epsilon array:'
2026           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2027           write (iout,'(/a)') 'One-body parameters:'
2028           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
2029               's||/s_|_^2','    chip  ','    alph  '
2030           write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2031                    sigii(i),chip(i),alp(i),i=1,ntyp)
2032         endif
2033        case default
2034         write(iout,*)"Wrong ipot"
2035 !   50 continue
2036       end select
2037       continue
2038       close (isidep)
2039 !-----------------------------------------------------------------------
2040 ! Calculate the "working" parameters of SC interactions.
2041
2042 !el from module energy - COMMON.INTERACT-------
2043       allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2044       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2045       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2046       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2047         dcavtub(ntyp1))
2048       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2049         tubetranene(ntyp1))
2050       aa_aq(:,:)=0.0D0
2051       bb_aq(:,:)=0.0D0
2052       aa_lip(:,:)=0.0D0
2053       bb_lip(:,:)=0.0D0
2054       chi(:,:)=0.0D0
2055       sigma(:,:)=0.0D0
2056       r0(:,:)=0.0D0
2057       acavtub(:)=0.0d0
2058       bcavtub(:)=0.0d0
2059       ccavtub(:)=0.0d0
2060       dcavtub(:)=0.0d0
2061       sc_aa_tube_par(:)=0.0d0
2062       sc_bb_tube_par(:)=0.0d0
2063
2064 !--------------------------------
2065
2066       do i=2,ntyp
2067         do j=1,i-1
2068           eps(i,j)=eps(j,i)
2069           epslip(i,j)=epslip(j,i)
2070         enddo
2071       enddo
2072       do i=1,ntyp
2073         do j=i,ntyp
2074           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2075           sigma(j,i)=sigma(i,j)
2076           rs0(i,j)=dwa16*sigma(i,j)
2077           rs0(j,i)=rs0(i,j)
2078         enddo
2079       enddo
2080       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2081        'Working parameters of the SC interactions:',&
2082        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
2083        '  chi1   ','   chi2   '
2084       do i=1,ntyp
2085         do j=i,ntyp
2086           epsij=eps(i,j)
2087           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2088             rrij=sigma(i,j)
2089           else
2090             rrij=rr0(i)+rr0(j)
2091           endif
2092           r0(i,j)=rrij
2093           r0(j,i)=rrij
2094           rrij=rrij**expon
2095           epsij=eps(i,j)
2096           sigeps=dsign(1.0D0,epsij)
2097           epsij=dabs(epsij)
2098           aa_aq(i,j)=epsij*rrij*rrij
2099           bb_aq(i,j)=-sigeps*epsij*rrij
2100           aa_aq(j,i)=aa_aq(i,j)
2101           bb_aq(j,i)=bb_aq(i,j)
2102           epsijlip=epslip(i,j)
2103           sigeps=dsign(1.0D0,epsijlip)
2104           epsijlip=dabs(epsijlip)
2105           aa_lip(i,j)=epsijlip*rrij*rrij
2106           bb_lip(i,j)=-sigeps*epsijlip*rrij
2107           aa_lip(j,i)=aa_lip(i,j)
2108           bb_lip(j,i)=bb_lip(i,j)
2109 !C          write(iout,*) aa_lip
2110           if (ipot.gt.2) then
2111             sigt1sq=sigma0(i)**2
2112             sigt2sq=sigma0(j)**2
2113             sigii1=sigii(i)
2114             sigii2=sigii(j)
2115             ratsig1=sigt2sq/sigt1sq
2116             ratsig2=1.0D0/ratsig1
2117             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2118             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2119             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2120           else
2121             rsum_max=sigma(i,j)
2122           endif
2123 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2124             sigmaii(i,j)=rsum_max
2125             sigmaii(j,i)=rsum_max
2126 !         else
2127 !           sigmaii(i,j)=r0(i,j)
2128 !           sigmaii(j,i)=r0(i,j)
2129 !         endif
2130 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2131           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2132             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2133             augm(i,j)=epsij*r_augm**(2*expon)
2134 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2135             augm(j,i)=augm(i,j)
2136           else
2137             augm(i,j)=0.0D0
2138             augm(j,i)=0.0D0
2139           endif
2140           if (lprint) then
2141             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2142             restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2143             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2144           endif
2145         enddo
2146       enddo
2147       write(iout,*) "tube param"
2148       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2149       ccavtubpep,dcavtubpep,tubetranenepep
2150       sigmapeptube=sigmapeptube**6
2151       sigeps=dsign(1.0D0,epspeptube)
2152       epspeptube=dabs(epspeptube)
2153       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2154       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2155       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2156       do i=1,ntyp
2157        read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2158       ccavtub(i),dcavtub(i),tubetranene(i)
2159        sigmasctube=sigmasctube**6
2160        sigeps=dsign(1.0D0,epssctube)
2161        epssctube=dabs(epssctube)
2162        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2163        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2164       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2165       enddo
2166
2167       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2168       bad(:,:)=0.0D0
2169
2170 #ifdef OLDSCP
2171 !
2172 ! Define the SC-p interaction constants (hard-coded; old style)
2173 !
2174       do i=1,ntyp
2175 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2176 ! helix formation)
2177 !       aad(i,1)=0.3D0*4.0D0**12
2178 ! Following line for constants currently implemented
2179 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2180         aad(i,1)=1.5D0*4.0D0**12
2181 !       aad(i,1)=0.17D0*5.6D0**12
2182         aad(i,2)=aad(i,1)
2183 ! "Soft" SC-p repulsion
2184         bad(i,1)=0.0D0
2185 ! Following line for constants currently implemented
2186 !       aad(i,1)=0.3D0*4.0D0**6
2187 ! "Hard" SC-p repulsion
2188         bad(i,1)=3.0D0*4.0D0**6
2189 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2190         bad(i,2)=bad(i,1)
2191 !       aad(i,1)=0.0D0
2192 !       aad(i,2)=0.0D0
2193 !       bad(i,1)=1228.8D0
2194 !       bad(i,2)=1228.8D0
2195       enddo
2196 #else
2197 !
2198 ! 8/9/01 Read the SC-p interaction constants from file
2199 !
2200       do i=1,ntyp
2201         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2202       enddo
2203       do i=1,ntyp
2204         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2205         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2206         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2207         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2208       enddo
2209 !      lprint=.true.
2210       if (lprint) then
2211         write (iout,*) "Parameters of SC-p interactions:"
2212         do i=1,ntyp
2213           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2214            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2215         enddo
2216       endif
2217 !      lprint=.false.
2218 #endif
2219 !
2220 ! Define the constants of the disulfide bridge
2221 !
2222       ebr=-5.50D0
2223 !
2224 ! Old arbitrary potential - commented out.
2225 !
2226 !      dbr= 4.20D0
2227 !      fbr= 3.30D0
2228 !
2229 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2230 ! energy surface of diethyl disulfide.
2231 ! A. Liwo and U. Kozlowska, 11/24/03
2232 !
2233       D0CM = 3.78d0
2234       AKCM = 15.1d0
2235       AKTH = 11.0d0
2236       AKCT = 12.0d0
2237       V1SS =-1.08d0
2238       V2SS = 7.61d0
2239       V3SS = 13.7d0
2240 !      akcm=0.0d0
2241 !      akth=0.0d0
2242 !      akct=0.0d0
2243 !      v1ss=0.0d0
2244 !      v2ss=0.0d0
2245 !      v3ss=0.0d0
2246
2247       if(me.eq.king) then
2248       write (iout,'(/a)') "Disulfide bridge parameters:"
2249       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2250       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2251       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2252       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2253         ' v3ss:',v3ss
2254       endif
2255       return
2256   111 write (iout,*) "Error reading bending energy parameters."
2257       goto 999
2258   112 write (iout,*) "Error reading rotamer energy parameters."
2259       goto 999
2260   113 write (iout,*) "Error reading torsional energy parameters."
2261       goto 999
2262   114 write (iout,*) "Error reading double torsional energy parameters."
2263       goto 999
2264   115 write (iout,*) &
2265         "Error reading cumulant (multibody energy) parameters."
2266       goto 999
2267   116 write (iout,*) "Error reading electrostatic energy parameters."
2268       goto 999
2269   117 write (iout,*) "Error reading side chain interaction parameters."
2270       goto 999
2271   118 write (iout,*) "Error reading SCp interaction parameters."
2272       goto 999
2273   119 write (iout,*) "Error reading SCCOR parameters"
2274   999 continue
2275 #ifdef MPI
2276       call MPI_Finalize(Ierror)
2277 #endif
2278       stop
2279       return
2280       end subroutine parmread
2281 #endif
2282 !-----------------------------------------------------------------------------
2283 ! printmat.f
2284 !-----------------------------------------------------------------------------
2285       subroutine printmat(ldim,m,n,iout,key,a)
2286
2287       integer :: n,ldim
2288       character(len=3),dimension(n) :: key
2289       real(kind=8),dimension(ldim,n) :: a
2290 !el local variables
2291       integer :: i,j,k,m,iout,nlim
2292
2293       do 1 i=1,n,8
2294       nlim=min0(i+7,n)
2295       write (iout,1000) (key(k),k=i,nlim)
2296       write (iout,1020)
2297  1000 format (/5x,8(6x,a3))
2298  1020 format (/80(1h-)/)
2299       do 2 j=1,n
2300       write (iout,1010) key(j),(a(j,k),k=i,nlim)
2301     2 continue
2302     1 continue
2303  1010 format (a3,2x,8(f9.4))
2304       return
2305       end subroutine printmat
2306 !-----------------------------------------------------------------------------
2307 ! readpdb.F
2308 !-----------------------------------------------------------------------------
2309       subroutine readpdb
2310 ! Read the PDB file and convert the peptide geometry into virtual-chain
2311 ! geometry.
2312       use geometry_data
2313       use energy_data, only: itype
2314       use control_data
2315       use compare_data
2316       use MPI_data
2317       use control, only: rescode
2318 !      implicit real*8 (a-h,o-z)
2319 !      include 'DIMENSIONS'
2320 !      include 'COMMON.LOCAL'
2321 !      include 'COMMON.VAR'
2322 !      include 'COMMON.CHAIN'
2323 !      include 'COMMON.INTERACT'
2324 !      include 'COMMON.IOUNITS'
2325 !      include 'COMMON.GEO'
2326 !      include 'COMMON.NAMES'
2327 !      include 'COMMON.CONTROL'
2328 !      include 'COMMON.DISTFIT'
2329 !      include 'COMMON.SETUP'
2330       integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2331 !        ishift_pdb
2332       logical :: lprn=.true.,fail
2333       real(kind=8),dimension(3) :: e1,e2,e3
2334       real(kind=8) :: dcj,efree_temp
2335       character(len=3) :: seq,res
2336       character(len=5) :: atom
2337       character(len=80) :: card
2338       real(kind=8),dimension(3,20) :: sccor
2339       integer :: kkk,lll,icha,kupa      !rescode,
2340       real(kind=8) :: cou
2341 !el local varables
2342       integer,dimension(2,maxres/3) :: hfrag_alloc
2343       integer,dimension(4,maxres/3) :: bfrag_alloc
2344       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2345
2346       efree_temp=0.0d0
2347       ibeg=1
2348       ishift1=0
2349       ishift=0
2350 !      write (2,*) "UNRES_PDB",unres_pdb
2351       ires=0
2352       ires_old=0
2353       nres=0
2354       iii=0
2355       lsecondary=.false.
2356       nhfrag=0
2357       nbfrag=0
2358 !-----------------------------
2359       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2360       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2361
2362       do i=1,100000
2363         read (ipdbin,'(a80)',end=10) card
2364 !       write (iout,'(a)') card
2365         if (card(:5).eq.'HELIX') then
2366           nhfrag=nhfrag+1
2367           lsecondary=.true.
2368           read(card(22:25),*) hfrag(1,nhfrag)
2369           read(card(34:37),*) hfrag(2,nhfrag)
2370         endif
2371         if (card(:5).eq.'SHEET') then
2372           nbfrag=nbfrag+1
2373           lsecondary=.true.
2374           read(card(24:26),*) bfrag(1,nbfrag)
2375           read(card(35:37),*) bfrag(2,nbfrag)
2376 !rc----------------------------------------
2377 !rc  to be corrected !!!
2378           bfrag(3,nbfrag)=bfrag(1,nbfrag)
2379           bfrag(4,nbfrag)=bfrag(2,nbfrag)
2380 !rc----------------------------------------
2381         endif
2382         if (card(:3).eq.'END') then
2383           goto 10
2384         else if (card(:3).eq.'TER') then
2385 ! End current chain
2386           ires_old=ires+2
2387           ishift1=ishift1+1
2388           itype(ires_old)=ntyp1
2389           itype(ires_old-1)=ntyp1
2390           ibeg=2
2391 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2392           if (unres_pdb) then
2393             do j=1,3
2394               dc(j,ires)=sccor(j,iii)
2395             enddo
2396           else
2397             call sccenter(ires,iii,sccor)
2398 !          iii=0
2399           endif
2400           iii=0
2401         endif
2402 ! Read free energy
2403         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2404 ! Fish out the ATOM cards.
2405         if (index(card(1:4),'ATOM').gt.0) then
2406           read (card(12:16),*) atom
2407 !          write (iout,*) "! ",atom," !",ires
2408 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2409           read (card(23:26),*) ires
2410           read (card(18:20),'(a3)') res
2411 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2412 !     &      " ires_old",ires_old
2413 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2414 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2415           if (ires-ishift+ishift1.ne.ires_old) then
2416 ! Calculate the CM of the preceding residue.
2417 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2418             if (ibeg.eq.0) then
2419 !              write (iout,*) "Calculating sidechain center iii",iii
2420               if (unres_pdb) then
2421                 do j=1,3
2422                   dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2423                 enddo
2424               else
2425                 call sccenter(ires_old,iii,sccor)
2426               endif
2427               iii=0
2428             endif
2429 ! Start new residue.
2430             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2431               ires=ires_old
2432               cycle
2433             else if (ibeg.eq.1) then
2434               write (iout,*) "BEG ires",ires
2435               ishift=ires-1
2436               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2437                 ishift=ishift-1
2438                 itype(1)=ntyp1
2439               endif
2440               ires=ires-ishift+ishift1
2441               ires_old=ires
2442 !              write (iout,*) "ishift",ishift," ires",ires,&
2443 !               " ires_old",ires_old
2444               ibeg=0
2445             else if (ibeg.eq.2) then
2446 ! Start a new chain
2447               ishift=-ires_old+ires-1 !!!!!
2448               ishift1=ishift1-1    !!!!!
2449 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2450               ires=ires-ishift+ishift1
2451               ires_old=ires
2452               ibeg=0
2453             else
2454               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2455               ires=ires-ishift+ishift1
2456               ires_old=ires
2457             endif
2458             if (res.eq.'ACE' .or. res.eq.'NHE') then
2459               itype(ires)=10
2460             else
2461               itype(ires)=rescode(ires,res,0)
2462             endif
2463           else
2464             ires=ires-ishift+ishift1
2465           endif
2466 !          write (iout,*) "ires_old",ires_old," ires",ires
2467           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2468 !            ishift1=ishift1+1
2469           endif
2470 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2471           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2472              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2473             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2474 !            write (iout,*) "backbone ",atom
2475 #ifdef DEBUG
2476             write (iout,'(2i3,2x,a,3f8.3)') &
2477             ires,itype(ires),res,(c(j,ires),j=1,3)
2478 #endif
2479             iii=iii+1
2480             do j=1,3
2481               sccor(j,iii)=c(j,ires)
2482             enddo
2483 !            write (*,*) card(23:27),ires,itype(ires)
2484           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2485                    atom.ne.'N' .and. atom.ne.'C' .and. &
2486                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2487                    atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2488 !            write (iout,*) "sidechain ",atom
2489             iii=iii+1
2490             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2491           endif
2492         endif
2493       enddo
2494    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2495       if (ires.eq.0) return
2496 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2497 ! system
2498       nres=ires
2499       do i=2,nres-1
2500 !        write (iout,*) i,itype(i)
2501 !        if (itype(i).eq.ntyp1) then
2502 !          write (iout,*) "dummy",i,itype(i)
2503 !          do j=1,3
2504 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2505 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
2506 !            dc(j,i)=c(j,i)
2507 !          enddo
2508 !        endif
2509         if (itype(i).eq.ntyp1) then
2510          if (itype(i+1).eq.ntyp1) then
2511 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2512 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
2513 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
2514            if (unres_pdb) then
2515 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2516 !            print *,i,'tu dochodze'
2517             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2518             if (fail) then
2519               e2(1)=0.0d0
2520               e2(2)=1.0d0
2521               e2(3)=0.0d0
2522             endif !fail
2523             print *,i,'a tu?'
2524             do j=1,3
2525              c(j,i)=c(j,i-1)-1.9d0*e2(j)
2526             enddo
2527            else   !unres_pdb
2528            do j=1,3
2529              dcj=(c(j,i-2)-c(j,i-3))/2.0
2530             if (dcj.eq.0) dcj=1.23591524223
2531              c(j,i)=c(j,i-1)+dcj
2532              c(j,nres+i)=c(j,i)
2533            enddo
2534           endif   !unres_pdb
2535          else     !itype(i+1).eq.ntyp1
2536           if (unres_pdb) then
2537 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2538             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2539             if (fail) then
2540               e2(1)=0.0d0
2541               e2(2)=1.0d0
2542               e2(3)=0.0d0
2543             endif
2544             do j=1,3
2545               c(j,i)=c(j,i+1)-1.9d0*e2(j)
2546             enddo
2547           else !unres_pdb
2548            do j=1,3
2549             dcj=(c(j,i+3)-c(j,i+2))/2.0
2550             if (dcj.eq.0) dcj=1.23591524223
2551             c(j,i)=c(j,i+1)-dcj
2552             c(j,nres+i)=c(j,i)
2553            enddo
2554           endif !unres_pdb
2555          endif !itype(i+1).eq.ntyp1
2556         endif  !itype.eq.ntyp1
2557
2558       enddo
2559 ! Calculate the CM of the last side chain.
2560       if (iii.gt.0)  then
2561       if (unres_pdb) then
2562         do j=1,3
2563           dc(j,ires)=sccor(j,iii)
2564         enddo
2565       else
2566         call sccenter(ires,iii,sccor)
2567       endif
2568       endif
2569 !      nres=ires
2570       nsup=nres
2571       nstart_sup=1
2572       if (itype(nres).ne.10) then
2573         nres=nres+1
2574         itype(nres)=ntyp1
2575         if (unres_pdb) then
2576 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2577           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2578           if (fail) then
2579             e2(1)=0.0d0
2580             e2(2)=1.0d0
2581             e2(3)=0.0d0
2582           endif
2583           do j=1,3
2584             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2585           enddo
2586         else
2587         do j=1,3
2588           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2589           c(j,nres)=c(j,nres-1)+dcj
2590           c(j,2*nres)=c(j,nres)
2591         enddo
2592         endif
2593       endif
2594 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2595 #ifdef WHAM_RUN
2596       if (nres.ne.nres0) then
2597         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2598                        " NRES0=",nres0
2599         stop "Error nres value in WHAM input"
2600       endif
2601 #endif
2602 !---------------------------------
2603 !el reallocate tables
2604 !      do i=1,maxres/3
2605 !       do j=1,2
2606 !         hfrag_alloc(j,i)=hfrag(j,i)
2607 !        enddo
2608 !       do j=1,4
2609 !         bfrag_alloc(j,i)=bfrag(j,i)
2610 !        enddo
2611 !      enddo
2612
2613 !      deallocate(hfrag)
2614 !      deallocate(bfrag)
2615 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
2616 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2617 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2618 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
2619
2620 !      do i=1,nhfrag
2621 !       do j=1,2
2622 !         hfrag(j,i)=hfrag_alloc(j,i)
2623 !        enddo
2624 !      enddo
2625 !      do i=1,nbfrag
2626 !       do j=1,4
2627 !         bfrag(j,i)=bfrag_alloc(j,i)
2628 !        enddo
2629 !      enddo
2630 !el end reallocate tables
2631 !---------------------------------
2632       do i=2,nres-1
2633         do j=1,3
2634           c(j,i+nres)=dc(j,i)
2635         enddo
2636       enddo
2637       do j=1,3
2638         c(j,nres+1)=c(j,1)
2639         c(j,2*nres)=c(j,nres)
2640       enddo
2641       if (itype(1).eq.ntyp1) then
2642         nsup=nsup-1
2643         nstart_sup=2
2644         if (unres_pdb) then
2645 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2646           call refsys(2,3,4,e1,e2,e3,fail)
2647           if (fail) then
2648             e2(1)=0.0d0
2649             e2(2)=1.0d0
2650             e2(3)=0.0d0
2651           endif
2652           do j=1,3
2653             c(j,1)=c(j,2)-1.9d0*e2(j)
2654           enddo
2655         else
2656         do j=1,3
2657           dcj=(c(j,4)-c(j,3))/2.0
2658           c(j,1)=c(j,2)-dcj
2659           c(j,nres+1)=c(j,1)
2660         enddo
2661         endif
2662       endif
2663 ! Copy the coordinates to reference coordinates
2664 !      do i=1,2*nres
2665 !        do j=1,3
2666 !          cref(j,i)=c(j,i)
2667 !        enddo
2668 !      enddo
2669 ! Calculate internal coordinates.
2670       if (lprn) then
2671       write (iout,'(/a)') &
2672         "Cartesian coordinates of the reference structure"
2673       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2674        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2675       do ires=1,nres
2676         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2677           restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2678           (c(j,ires+nres),j=1,3)
2679       enddo
2680       endif
2681 ! znamy już nres wiec mozna alokowac tablice
2682 ! Calculate internal coordinates.
2683       if(me.eq.king.or..not.out1file)then
2684        write (iout,'(a)') &
2685          "Backbone and SC coordinates as read from the PDB"
2686        do ires=1,nres
2687         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2688           ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2689           (c(j,nres+ires),j=1,3)
2690        enddo
2691       endif
2692
2693       if(.not.allocated(vbld)) then
2694        allocate(vbld(2*nres))
2695        do i=1,2*nres
2696          vbld(i)=0.d0
2697        enddo
2698       endif
2699       if(.not.allocated(vbld_inv)) then
2700        allocate(vbld_inv(2*nres))
2701        do i=1,2*nres
2702          vbld_inv(i)=0.d0
2703        enddo
2704       endif
2705 !!!el
2706       if(.not.allocated(theta)) then
2707         allocate(theta(nres+2))
2708         theta(:)=0.0d0
2709       endif
2710
2711       if(.not.allocated(phi)) allocate(phi(nres+2))
2712       if(.not.allocated(alph)) allocate(alph(nres+2))
2713       if(.not.allocated(omeg)) allocate(omeg(nres+2))
2714       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2715       if(.not.allocated(phiref)) allocate(phiref(nres+2))
2716       if(.not.allocated(costtab)) allocate(costtab(nres))
2717       if(.not.allocated(sinttab)) allocate(sinttab(nres))
2718       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2719       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2720       if(.not.allocated(xxref)) allocate(xxref(nres))
2721       if(.not.allocated(yyref)) allocate(yyref(nres))
2722       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2723       if(.not.allocated(dc_norm)) then
2724 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2725         allocate(dc_norm(3,0:2*nres+2))
2726         dc_norm(:,:)=0.d0
2727       endif
2728
2729       call int_from_cart(.true.,.false.)
2730       call sc_loc_geom(.false.)
2731       do i=1,nres
2732         thetaref(i)=theta(i)
2733         phiref(i)=phi(i)
2734       enddo
2735 !      do i=1,2*nres
2736 !        vbld_inv(i)=0.d0
2737 !        vbld(i)=0.d0
2738 !      enddo
2739
2740       do i=1,nres-1
2741         do j=1,3
2742           dc(j,i)=c(j,i+1)-c(j,i)
2743           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2744         enddo
2745       enddo
2746       do i=2,nres-1
2747         do j=1,3
2748           dc(j,i+nres)=c(j,i+nres)-c(j,i)
2749           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2750         enddo
2751 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2752 !        vbld_inv(i+nres)
2753       enddo
2754 !      call chainbuild
2755 ! Copy the coordinates to reference coordinates
2756 ! Splits to single chain if occurs
2757
2758 !      do i=1,2*nres
2759 !        do j=1,3
2760 !          cref(j,i,cou)=c(j,i)
2761 !        enddo
2762 !      enddo
2763 !
2764       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2765       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2766       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2767 !-----------------------------
2768       kkk=1
2769       lll=0
2770       cou=1
2771       do i=1,nres
2772       lll=lll+1
2773 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2774       if (i.gt.1) then
2775       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2776       chain_length=lll-1
2777       kkk=kkk+1
2778 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2779       lll=1
2780       endif
2781       endif
2782         do j=1,3
2783           cref(j,i,cou)=c(j,i)
2784           cref(j,i+nres,cou)=c(j,i+nres)
2785           if (i.le.nres) then
2786           chain_rep(j,lll,kkk)=c(j,i)
2787           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2788           endif
2789          enddo
2790       enddo
2791       write (iout,*) chain_length
2792       if (chain_length.eq.0) chain_length=nres
2793       do j=1,3
2794       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2795       chain_rep(j,chain_length+nres,symetr) &
2796       =chain_rep(j,chain_length+nres,1)
2797       enddo
2798 ! diagnostic
2799 !       write (iout,*) "spraw lancuchy",chain_length,symetr
2800 !       do i=1,4
2801 !         do kkk=1,chain_length
2802 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2803 !         enddo
2804 !        enddo
2805 ! enddiagnostic
2806 ! makes copy of chains
2807         write (iout,*) "symetr", symetr
2808       do j=1,3
2809       dc(j,0)=c(j,1)
2810       enddo
2811
2812       if (symetr.gt.1) then
2813        call permut(symetr)
2814        nperm=1
2815        do i=1,symetr
2816        nperm=nperm*i
2817        enddo
2818        do i=1,nperm
2819        write(iout,*) (tabperm(i,kkk),kkk=1,4)
2820        enddo
2821        do i=1,nperm
2822         cou=0
2823         do kkk=1,symetr
2824          icha=tabperm(i,kkk)
2825 !         write (iout,*) i,icha
2826          do lll=1,chain_length
2827           cou=cou+1
2828            if (cou.le.nres) then
2829            do j=1,3
2830             kupa=mod(lll,chain_length)
2831             iprzes=(kkk-1)*chain_length+lll
2832             if (kupa.eq.0) kupa=chain_length
2833 !            write (iout,*) "kupa", kupa
2834             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2835             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2836           enddo
2837           endif
2838          enddo
2839         enddo
2840        enddo
2841        endif
2842 !-koniec robienia kopii
2843 ! diag
2844       do kkk=1,nperm
2845       write (iout,*) "nowa struktura", nperm
2846       do i=1,nres
2847         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2848       cref(2,i,kkk),&
2849       cref(3,i,kkk),cref(1,nres+i,kkk),&
2850       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2851       enddo
2852   100 format (//'              alpha-carbon coordinates       ',&
2853                 '     centroid coordinates'/ &
2854                 '       ', 6X,'X',11X,'Y',11X,'Z', &
2855                                 10X,'X',11X,'Y',11X,'Z')
2856   110 format (a,'(',i3,')',6f12.5)
2857
2858       enddo
2859 !c enddiag
2860       do j=1,nbfrag
2861         do i=1,4
2862          bfrag(i,j)=bfrag(i,j)-ishift
2863         enddo
2864       enddo
2865
2866       do j=1,nhfrag
2867         do i=1,2
2868          hfrag(i,j)=hfrag(i,j)-ishift
2869         enddo
2870       enddo
2871       ishift_pdb=ishift
2872
2873       return
2874       end subroutine readpdb
2875 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2876 !-----------------------------------------------------------------------------
2877 ! readrtns_CSA.F
2878 !-----------------------------------------------------------------------------
2879       subroutine read_control
2880 !
2881 ! Read contorl data
2882 !
2883 !      use geometry_data
2884       use comm_machsw
2885       use energy_data
2886       use control_data
2887       use compare_data
2888       use MCM_data
2889       use map_data
2890       use csa_data
2891       use MD_data
2892       use MPI_data
2893       use random, only: random_init
2894 !      implicit real*8 (a-h,o-z)
2895 !      include 'DIMENSIONS'
2896 #ifdef MP
2897       use prng, only:prng_restart
2898       include 'mpif.h'
2899       logical :: OKRandom!, prng_restart
2900       real(kind=8) :: r1
2901 #endif
2902 !      include 'COMMON.IOUNITS'
2903 !      include 'COMMON.TIME1'
2904 !      include 'COMMON.THREAD'
2905 !      include 'COMMON.SBRIDGE'
2906 !      include 'COMMON.CONTROL'
2907 !      include 'COMMON.MCM'
2908 !      include 'COMMON.MAP'
2909 !      include 'COMMON.HEADER'
2910 !      include 'COMMON.CSA'
2911 !      include 'COMMON.CHAIN'
2912 !      include 'COMMON.MUCA'
2913 !      include 'COMMON.MD'
2914 !      include 'COMMON.FFIELD'
2915 !      include 'COMMON.INTERACT'
2916 !      include 'COMMON.SETUP'
2917 !el      integer :: KDIAG,ICORFL,IXDR
2918 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2919       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2920         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
2921 !      character(len=80) :: ucase
2922       character(len=640) :: controlcard
2923
2924       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2925       integer i
2926
2927       nglob_csa=0
2928       eglob_csa=1d99
2929       nmin_csa=0
2930       read (INP,'(a)') titel
2931       call card_concat(controlcard,.true.)
2932 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2933 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2934       call reada(controlcard,'SEED',seed,0.0D0)
2935       call random_init(seed)
2936 ! Set up the time limit (caution! The time must be input in minutes!)
2937       read_cart=index(controlcard,'READ_CART').gt.0
2938       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2939       call readi(controlcard,'SYM',symetr,1)
2940       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2941       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2942       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2943       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2944       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2945       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2946       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2947       call reada(controlcard,'DRMS',drms,0.1D0)
2948       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2949        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2950        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2951        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2952        write (iout,'(a,f10.1)')'DRMS    = ',drms
2953        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2954        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2955       endif
2956       call readi(controlcard,'NZ_START',nz_start,0)
2957       call readi(controlcard,'NZ_END',nz_end,0)
2958       call readi(controlcard,'IZ_SC',iz_sc,0)
2959       timlim=60.0D0*timlim
2960       safety = 60.0d0*safety
2961       timem=timlim
2962       modecalc=0
2963       call reada(controlcard,"T_BATH",t_bath,300.0d0)
2964 !C SHIELD keyword sets if the shielding effect of side-chains is used
2965 !C 0 denots no shielding is used all peptide are equally despite the
2966 !C solvent accesible area
2967 !C 1 the newly introduced function
2968 !C 2 reseved for further possible developement
2969       call readi(controlcard,'SHIELD',shield_mode,0)
2970 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2971         write(iout,*) "shield_mode",shield_mode
2972 !C  Varibles set size of box
2973       call reada(controlcard,'BOXX',boxxsize,100.0d0)
2974       call reada(controlcard,'BOXY',boxysize,100.0d0)
2975       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2976       call readi(controlcard,'TUBEMOD',tubemode,0)
2977       write (iout,*) TUBEmode,"TUBEMODE"
2978       if (TUBEmode.gt.0) then
2979        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
2980        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
2981        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
2982        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
2983        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
2984        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
2985        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
2986        buftubebot=bordtubebot+tubebufthick
2987        buftubetop=bordtubetop-tubebufthick
2988       endif
2989
2990 ! CUTOFFF ON ELECTROSTATICS
2991       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2992       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2993       write(iout,*) "R_CUT_ELE=",r_cut_ele
2994 ! Lipidic parameters
2995       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
2996       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
2997       if (lipthick.gt.0.0d0) then
2998        bordliptop=(boxzsize+lipthick)/2.0
2999        bordlipbot=bordliptop-lipthick
3000       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3001       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3002       buflipbot=bordlipbot+lipbufthick
3003       bufliptop=bordliptop-lipbufthick
3004       if ((lipbufthick*2.0d0).gt.lipthick) &
3005        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3006       endif !lipthick.gt.0
3007       write(iout,*) "bordliptop=",bordliptop
3008       write(iout,*) "bordlipbot=",bordlipbot
3009       write(iout,*) "bufliptop=",bufliptop
3010       write(iout,*) "buflipbot=",buflipbot
3011       write (iout,*) "SHIELD MODE",shield_mode
3012
3013 !C-------------------------
3014       minim=(index(controlcard,'MINIMIZE').gt.0)
3015       dccart=(index(controlcard,'CART').gt.0)
3016       overlapsc=(index(controlcard,'OVERLAP').gt.0)
3017       overlapsc=.not.overlapsc
3018       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3019       searchsc=.not.searchsc
3020       sideadd=(index(controlcard,'SIDEADD').gt.0)
3021       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3022       outpdb=(index(controlcard,'PDBOUT').gt.0)
3023       outmol2=(index(controlcard,'MOL2OUT').gt.0)
3024       pdbref=(index(controlcard,'PDBREF').gt.0)
3025       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3026       indpdb=index(controlcard,'PDBSTART')
3027       extconf=(index(controlcard,'EXTCONF').gt.0)
3028       call readi(controlcard,'IPRINT',iprint,0)
3029       call readi(controlcard,'MAXGEN',maxgen,10000)
3030       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3031       call readi(controlcard,"KDIAG",kdiag,0)
3032       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3033       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3034        write (iout,*) "RESCALE_MODE",rescale_mode
3035       split_ene=index(controlcard,'SPLIT_ENE').gt.0
3036       if (index(controlcard,'REGULAR').gt.0.0D0) then
3037         call reada(controlcard,'WEIDIS',weidis,0.1D0)
3038         modecalc=1
3039         refstr=.true.
3040       endif
3041       if (index(controlcard,'CHECKGRAD').gt.0) then
3042         modecalc=5
3043         if (index(controlcard,'CART').gt.0) then
3044           icheckgrad=1
3045         elseif (index(controlcard,'CARINT').gt.0) then
3046           icheckgrad=2
3047         else
3048           icheckgrad=3
3049         endif
3050       elseif (index(controlcard,'THREAD').gt.0) then
3051         modecalc=2
3052         call readi(controlcard,'THREAD',nthread,0)
3053         if (nthread.gt.0) then
3054           call reada(controlcard,'WEIDIS',weidis,0.1D0)
3055         else
3056           if (fg_rank.eq.0) &
3057           write (iout,'(a)')'A number has to follow the THREAD keyword.'
3058           stop 'Error termination in Read_Control.'
3059         endif
3060       else if (index(controlcard,'MCMA').gt.0) then
3061         modecalc=3
3062       else if (index(controlcard,'MCEE').gt.0) then
3063         modecalc=6
3064       else if (index(controlcard,'MULTCONF').gt.0) then
3065         modecalc=4
3066       else if (index(controlcard,'MAP').gt.0) then
3067         modecalc=7
3068         call readi(controlcard,'MAP',nmap,0)
3069       else if (index(controlcard,'CSA').gt.0) then
3070         modecalc=8
3071 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3072 !rc
3073 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3074 !rc
3075 !rc        modecalc=9
3076 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3077 !fmc        modecalc=10
3078       else if (index(controlcard,'SOFTREG').gt.0) then
3079         modecalc=11
3080       else if (index(controlcard,'CHECK_BOND').gt.0) then
3081         modecalc=-1
3082       else if (index(controlcard,'TEST').gt.0) then
3083         modecalc=-2
3084       else if (index(controlcard,'MD').gt.0) then
3085         modecalc=12
3086       else if (index(controlcard,'RE ').gt.0) then
3087         modecalc=14
3088       endif
3089
3090       lmuca=index(controlcard,'MUCA').gt.0
3091       call readi(controlcard,'MUCADYN',mucadyn,0)
3092       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3093       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3094        then
3095        write (iout,*) 'MUCADYN=',mucadyn
3096        write (iout,*) 'MUCASMOOTH=',muca_smooth
3097       endif
3098
3099       iscode=index(controlcard,'ONE_LETTER')
3100       indphi=index(controlcard,'PHI')
3101       indback=index(controlcard,'BACK')
3102       iranconf=index(controlcard,'RAND_CONF')
3103       i2ndstr=index(controlcard,'USE_SEC_PRED')
3104       gradout=index(controlcard,'GRADOUT').gt.0
3105       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3106       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3107       if (me.eq.king .or. .not.out1file ) &
3108         write (iout,*) "DISTCHAINMAX",distchainmax
3109
3110       if(me.eq.king.or..not.out1file) &
3111        write (iout,'(2a)') diagmeth(kdiag),&
3112         ' routine used to diagonalize matrices.'
3113       if (shield_mode.gt.0) then
3114       pi=3.141592d0
3115 !C VSolvSphere the volume of solving sphere
3116 !C      print *,pi,"pi"
3117 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3118 !C there will be no distinction between proline peptide group and normal peptide
3119 !C group in case of shielding parameters
3120       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3121       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3122       write (iout,*) VSolvSphere,VSolvSphere_div
3123 !C long axis of side chain
3124       do i=1,ntyp
3125       long_r_sidechain(i)=vbldsc0(1,i)
3126       short_r_sidechain(i)=sigma0(i)
3127       write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3128          sigma0(i)
3129       enddo
3130       buff_shield=1.0d0
3131       endif
3132       return
3133       end subroutine read_control
3134 !-----------------------------------------------------------------------------
3135       subroutine read_REMDpar
3136 !
3137 ! Read REMD settings
3138 !
3139 !       use control
3140 !       use energy
3141 !       use geometry
3142       use REMD_data
3143       use MPI_data
3144       use control_data, only:out1file
3145 !      implicit real*8 (a-h,o-z)
3146 !      include 'DIMENSIONS'
3147 !      include 'COMMON.IOUNITS'
3148 !      include 'COMMON.TIME1'
3149 !      include 'COMMON.MD'
3150       use MD_data
3151 !el #ifndef LANG0
3152 !el      include 'COMMON.LANGEVIN'
3153 !el #else
3154 !el      include 'COMMON.LANGEVIN.lang0'
3155 !el #endif
3156 !      include 'COMMON.INTERACT'
3157 !      include 'COMMON.NAMES'
3158 !      include 'COMMON.GEO'
3159 !      include 'COMMON.REMD'
3160 !      include 'COMMON.CONTROL'
3161 !      include 'COMMON.SETUP'
3162 !      character(len=80) :: ucase
3163       character(len=320) :: controlcard
3164       character(len=3200) :: controlcard1
3165       integer :: iremd_m_total
3166 !el local variables
3167       integer :: i
3168 !     real(kind=8) :: var,ene
3169
3170       if(me.eq.king.or..not.out1file) &
3171        write (iout,*) "REMD setup"
3172
3173       call card_concat(controlcard,.true.)
3174       call readi(controlcard,"NREP",nrep,3)
3175       call readi(controlcard,"NSTEX",nstex,1000)
3176       call reada(controlcard,"RETMIN",retmin,10.0d0)
3177       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3178       mremdsync=(index(controlcard,'SYNC').gt.0)
3179       call readi(controlcard,"NSYN",i_sync_step,100)
3180       restart1file=(index(controlcard,'REST1FILE').gt.0)
3181       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3182       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3183       if(max_cache_traj_use.gt.max_cache_traj) &
3184                 max_cache_traj_use=max_cache_traj
3185       if(me.eq.king.or..not.out1file) then
3186 !d       if (traj1file) then
3187 !rc caching is in testing - NTWX is not ignored
3188 !d        write (iout,*) "NTWX value is ignored"
3189 !d        write (iout,*) "  trajectory is stored to one file by master"
3190 !d        write (iout,*) "  before exchange at NSTEX intervals"
3191 !d       endif
3192        write (iout,*) "NREP= ",nrep
3193        write (iout,*) "NSTEX= ",nstex
3194        write (iout,*) "SYNC= ",mremdsync
3195        write (iout,*) "NSYN= ",i_sync_step
3196        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3197       endif
3198       remd_tlist=.false.
3199       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3200       if (index(controlcard,'TLIST').gt.0) then
3201          remd_tlist=.true.
3202          call card_concat(controlcard1,.true.)
3203          read(controlcard1,*) (remd_t(i),i=1,nrep)
3204          if(me.eq.king.or..not.out1file) &
3205           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3206       endif
3207       remd_mlist=.false.
3208       if (index(controlcard,'MLIST').gt.0) then
3209          remd_mlist=.true.
3210          call card_concat(controlcard1,.true.)
3211          read(controlcard1,*) (remd_m(i),i=1,nrep)
3212          if(me.eq.king.or..not.out1file) then
3213           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3214           iremd_m_total=0
3215           do i=1,nrep
3216            iremd_m_total=iremd_m_total+remd_m(i)
3217           enddo
3218           write (iout,*) 'Total number of replicas ',iremd_m_total
3219          endif
3220       endif
3221       if(me.eq.king.or..not.out1file) &
3222        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3223       return
3224       end subroutine read_REMDpar
3225 !-----------------------------------------------------------------------------
3226       subroutine read_MDpar
3227 !
3228 ! Read MD settings
3229 !
3230       use control_data, only: r_cut,rlamb,out1file
3231       use energy_data
3232       use geometry_data, only: pi
3233       use MPI_data
3234 !      implicit real*8 (a-h,o-z)
3235 !      include 'DIMENSIONS'
3236 !      include 'COMMON.IOUNITS'
3237 !      include 'COMMON.TIME1'
3238 !      include 'COMMON.MD'
3239       use MD_data
3240 !el #ifndef LANG0
3241 !el      include 'COMMON.LANGEVIN'
3242 !el #else
3243 !el      include 'COMMON.LANGEVIN.lang0'
3244 !el #endif
3245 !      include 'COMMON.INTERACT'
3246 !      include 'COMMON.NAMES'
3247 !      include 'COMMON.GEO'
3248 !      include 'COMMON.SETUP'
3249 !      include 'COMMON.CONTROL'
3250 !      include 'COMMON.SPLITELE'
3251 !      character(len=80) :: ucase
3252       character(len=320) :: controlcard
3253 !el local variables
3254       integer :: i
3255       real(kind=8) :: eta
3256
3257       call card_concat(controlcard,.true.)
3258       call readi(controlcard,"NSTEP",n_timestep,1000000)
3259       call readi(controlcard,"NTWE",ntwe,100)
3260       call readi(controlcard,"NTWX",ntwx,1000)
3261       call reada(controlcard,"DT",d_time,1.0d-1)
3262       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3263       call reada(controlcard,"DAMAX",damax,1.0d1)
3264       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3265       call readi(controlcard,"LANG",lang,0)
3266       RESPA = index(controlcard,"RESPA") .gt. 0
3267       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3268       ntime_split0=ntime_split
3269       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3270       ntime_split0=ntime_split
3271       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3272       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3273       rest = index(controlcard,"REST").gt.0
3274       tbf = index(controlcard,"TBF").gt.0
3275       usampl = index(controlcard,"USAMPL").gt.0
3276       mdpdb = index(controlcard,"MDPDB").gt.0
3277       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3278       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3279       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3280       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3281       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3282       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3283       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3284       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3285       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3286       large = index(controlcard,"LARGE").gt.0
3287       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3288       rattle = index(controlcard,"RATTLE").gt.0
3289 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3290       nset=0
3291       if(usampl) then
3292         call read_fragments
3293       endif
3294
3295       if(me.eq.king.or..not.out1file) then
3296        write (iout,*)
3297        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3298        write (iout,*)
3299        write (iout,'(a)') "The units are:"
3300        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3301        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3302         " acceleration: angstrom/(48.9 fs)**2"
3303        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3304        write (iout,*)
3305        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3306        write (iout,'(a60,f10.5,a)') &
3307         "Initial time step of numerical integration:",d_time,&
3308         " natural units"
3309        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3310        if (RESPA) then
3311         write (iout,'(2a,i4,a)') &
3312           "A-MTS algorithm used; initial time step for fast-varying",&
3313           " short-range forces split into",ntime_split," steps."
3314         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3315          r_cut," lambda",rlamb
3316        endif
3317        write (iout,'(2a,f10.5)') &
3318         "Maximum acceleration threshold to reduce the time step",&
3319         "/increase split number:",damax
3320        write (iout,'(2a,f10.5)') &
3321         "Maximum predicted energy drift to reduce the timestep",&
3322         "/increase split number:",edriftmax
3323        write (iout,'(a60,f10.5)') &
3324        "Maximum velocity threshold to reduce velocities:",dvmax
3325        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3326        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3327        if (rattle) write (iout,'(a60)') &
3328         "Rattle algorithm used to constrain the virtual bonds"
3329       endif
3330       reset_fricmat=1000
3331       if (lang.gt.0) then
3332         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3333         call reada(controlcard,"RWAT",rwat,1.4d0)
3334         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3335         surfarea=index(controlcard,"SURFAREA").gt.0
3336         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3337         if(me.eq.king.or..not.out1file)then
3338          write (iout,'(/a,$)') "Langevin dynamics calculation"
3339          if (lang.eq.1) then
3340           write (iout,'(a/)') &
3341             " with direct integration of Langevin equations"
3342          else if (lang.eq.2) then
3343           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3344          else if (lang.eq.3) then
3345           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3346          else if (lang.eq.4) then
3347           write (iout,'(a/)') " in overdamped mode"
3348          else
3349           write (iout,'(//a,i5)') &
3350             "=========== ERROR: Unknown Langevin dynamics mode:",lang
3351           stop
3352          endif
3353          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3354          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3355          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3356          write (iout,'(a60,f10.5)') &
3357          "Scaling factor of the friction forces:",scal_fric
3358          if (surfarea) write (iout,'(2a,i10,a)') &
3359            "Friction coefficients will be scaled by solvent-accessible",&
3360            " surface area every",reset_fricmat," steps."
3361         endif
3362 ! Calculate friction coefficients and bounds of stochastic forces
3363         eta=6*pi*cPoise*etawat
3364         if(me.eq.king.or..not.out1file) &
3365          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3366           eta
3367         gamp=scal_fric*(pstok+rwat)*eta
3368         stdfp=dsqrt(2*Rb*t_bath/d_time)
3369         allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3370         do i=1,ntyp
3371           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3372           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3373         enddo
3374         if(me.eq.king.or..not.out1file)then
3375          write (iout,'(/2a/)') &
3376          "Radii of site types and friction coefficients and std's of",&
3377          " stochastic forces of fully exposed sites"
3378          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3379          do i=1,ntyp
3380           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3381            gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3382          enddo
3383         endif
3384       else if (tbf) then
3385         if(me.eq.king.or..not.out1file)then
3386          write (iout,'(a)') "Berendsen bath calculation"
3387          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3388          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3389          if (reset_moment) &
3390          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3391          count_reset_moment," steps"
3392          if (reset_vel) &
3393           write (iout,'(a,i10,a)') &
3394           "Velocities will be reset at random every",count_reset_vel,&
3395          " steps"
3396         endif
3397       else
3398         if(me.eq.king.or..not.out1file) &
3399          write (iout,'(a31)') "Microcanonical mode calculation"
3400       endif
3401       if(me.eq.king.or..not.out1file)then
3402        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3403        if (usampl) then
3404           write(iout,*) "MD running with constraints."
3405           write(iout,*) "Equilibration time ", eq_time, " mtus."
3406           write(iout,*) "Constraining ", nfrag," fragments."
3407           write(iout,*) "Length of each fragment, weight and q0:"
3408           do iset=1,nset
3409            write (iout,*) "Set of restraints #",iset
3410            do i=1,nfrag
3411               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3412                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3413            enddo
3414            write(iout,*) "constraints between ", npair, "fragments."
3415            write(iout,*) "constraint pairs, weights and q0:"
3416            do i=1,npair
3417             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3418                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3419            enddo
3420            write(iout,*) "angle constraints within ", nfrag_back,&
3421             "backbone fragments."
3422            write(iout,*) "fragment, weights:"
3423            do i=1,nfrag_back
3424             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3425                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3426                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3427            enddo
3428           enddo
3429         iset=mod(kolor,nset)+1
3430        endif
3431       endif
3432       if(me.eq.king.or..not.out1file) &
3433        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3434       return
3435       end subroutine read_MDpar
3436 !-----------------------------------------------------------------------------
3437       subroutine map_read
3438
3439       use map_data
3440 !      implicit real*8 (a-h,o-z)
3441 !      include 'DIMENSIONS'
3442 !      include 'COMMON.MAP'
3443 !      include 'COMMON.IOUNITS'
3444       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3445       character(len=80) :: mapcard      !,ucase
3446 !el local variables
3447       integer :: imap
3448 !     real(kind=8) :: var,ene
3449
3450       do imap=1,nmap
3451         read (inp,'(a)') mapcard
3452         mapcard=ucase(mapcard)
3453         if (index(mapcard,'PHI').gt.0) then
3454           kang(imap)=1
3455         else if (index(mapcard,'THE').gt.0) then
3456           kang(imap)=2
3457         else if (index(mapcard,'ALP').gt.0) then
3458           kang(imap)=3
3459         else if (index(mapcard,'OME').gt.0) then
3460           kang(imap)=4
3461         else
3462           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3463           stop 'Error - illegal variable spec in MAP card.'
3464         endif
3465         call readi (mapcard,'RES1',res1(imap),0)
3466         call readi (mapcard,'RES2',res2(imap),0)
3467         if (res1(imap).eq.0) then
3468           res1(imap)=res2(imap)
3469         else if (res2(imap).eq.0) then
3470           res2(imap)=res1(imap)
3471         endif
3472         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3473           write (iout,'(a)') &
3474           'Error - illegal definition of variable group in MAP.'
3475           stop 'Error - illegal definition of variable group in MAP.'
3476         endif
3477         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3478         call reada(mapcard,'TO',ang_to(imap),0.0D0)
3479         call readi(mapcard,'NSTEP',nstep(imap),0)
3480         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3481           write (iout,'(a)') &
3482            'Illegal boundary and/or step size specification in MAP.'
3483           stop 'Illegal boundary and/or step size specification in MAP.'
3484         endif
3485       enddo ! imap
3486       return
3487       end subroutine map_read
3488 !-----------------------------------------------------------------------------
3489       subroutine csaread
3490
3491       use control_data, only: vdisulf
3492       use csa_data
3493 !      implicit real*8 (a-h,o-z)
3494 !      include 'DIMENSIONS'
3495 !      include 'COMMON.IOUNITS'
3496 !      include 'COMMON.GEO'
3497 !      include 'COMMON.CSA'
3498 !      include 'COMMON.BANK'
3499 !      include 'COMMON.CONTROL'
3500 !      character(len=80) :: ucase
3501       character(len=620) :: mcmcard
3502 !el local variables
3503 !     integer :: ntf,ik,iw_pdb
3504 !     real(kind=8) :: var,ene
3505
3506       call card_concat(mcmcard,.true.)
3507
3508       call readi(mcmcard,'NCONF',nconf,50)
3509       call readi(mcmcard,'NADD',nadd,0)
3510       call readi(mcmcard,'JSTART',jstart,1)
3511       call readi(mcmcard,'JEND',jend,1)
3512       call readi(mcmcard,'NSTMAX',nstmax,500000)
3513       call readi(mcmcard,'N0',n0,1)
3514       call readi(mcmcard,'N1',n1,6)
3515       call readi(mcmcard,'N2',n2,4)
3516       call readi(mcmcard,'N3',n3,0)
3517       call readi(mcmcard,'N4',n4,0)
3518       call readi(mcmcard,'N5',n5,0)
3519       call readi(mcmcard,'N6',n6,10)
3520       call readi(mcmcard,'N7',n7,0)
3521       call readi(mcmcard,'N8',n8,0)
3522       call readi(mcmcard,'N9',n9,0)
3523       call readi(mcmcard,'N14',n14,0)
3524       call readi(mcmcard,'N15',n15,0)
3525       call readi(mcmcard,'N16',n16,0)
3526       call readi(mcmcard,'N17',n17,0)
3527       call readi(mcmcard,'N18',n18,0)
3528
3529       vdisulf=(index(mcmcard,'DYNSS').gt.0)
3530
3531       call readi(mcmcard,'NDIFF',ndiff,2)
3532       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3533       call readi(mcmcard,'IS1',is1,1)
3534       call readi(mcmcard,'IS2',is2,8)
3535       call readi(mcmcard,'NRAN0',nran0,4)
3536       call readi(mcmcard,'NRAN1',nran1,2)
3537       call readi(mcmcard,'IRR',irr,1)
3538       call readi(mcmcard,'NSEED',nseed,20)
3539       call readi(mcmcard,'NTOTAL',ntotal,10000)
3540       call reada(mcmcard,'CUT1',cut1,2.0d0)
3541       call reada(mcmcard,'CUT2',cut2,5.0d0)
3542       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3543       call readi(mcmcard,'ICMAX',icmax,3)
3544       call readi(mcmcard,'IRESTART',irestart,0)
3545 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
3546       ntbankm=0
3547 !!bankt
3548       call reada(mcmcard,'DELE',dele,20.0d0)
3549       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3550       call readi(mcmcard,'IREF',iref,0)
3551       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3552       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3553       call readi(mcmcard,'NCONF_IN',nconf_in,0)
3554       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3555       write (iout,*) "NCONF_IN",nconf_in
3556       return
3557       end subroutine csaread
3558 !-----------------------------------------------------------------------------
3559       subroutine mcmread
3560
3561       use mcm_data
3562       use control_data, only: MaxMoveType
3563       use MD_data
3564       use minim_data
3565 !      implicit real*8 (a-h,o-z)
3566 !      include 'DIMENSIONS'
3567 !      include 'COMMON.MCM'
3568 !      include 'COMMON.MCE'
3569 !      include 'COMMON.IOUNITS'
3570 !      character(len=80) :: ucase
3571       character(len=320) :: mcmcard
3572 !el local variables
3573       integer :: i
3574 !     real(kind=8) :: var,ene
3575
3576       call card_concat(mcmcard,.true.)
3577       call readi(mcmcard,'MAXACC',maxacc,100)
3578       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3579       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3580       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3581       call readi(mcmcard,'MAXREPM',maxrepm,200)
3582       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3583       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3584       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3585       call reada(mcmcard,'E_UP',e_up,5.0D0)
3586       call reada(mcmcard,'DELTE',delte,0.1D0)
3587       call readi(mcmcard,'NSWEEP',nsweep,5)
3588       call readi(mcmcard,'NSTEPH',nsteph,0)
3589       call readi(mcmcard,'NSTEPC',nstepc,0)
3590       call reada(mcmcard,'TMIN',tmin,298.0D0)
3591       call reada(mcmcard,'TMAX',tmax,298.0D0)
3592       call readi(mcmcard,'NWINDOW',nwindow,0)
3593       call readi(mcmcard,'PRINT_MC',print_mc,0)
3594       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3595       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3596       ent_read=(index(mcmcard,'ENT_READ').gt.0)
3597       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3598       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3599       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3600       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3601       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3602       if (nwindow.gt.0) then
3603         allocate(winstart(nwindow))     !!el (maxres)
3604         allocate(winend(nwindow))       !!el
3605         allocate(winlen(nwindow))       !!el
3606         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3607         do i=1,nwindow
3608           winlen(i)=winend(i)-winstart(i)+1
3609         enddo
3610       endif
3611       if (tmax.lt.tmin) tmax=tmin
3612       if (tmax.eq.tmin) then
3613         nstepc=0
3614         nsteph=0
3615       endif
3616       if (nstepc.gt.0 .and. nsteph.gt.0) then
3617         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3618         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3619       endif
3620       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3621 ! Probabilities of different move types
3622       sumpro_type(0)=0.0D0
3623       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3624       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3625       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3626       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
3627       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3628       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3629       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3630       do i=1,MaxMoveType
3631         print *,'i',i,' sumprotype',sumpro_type(i)
3632         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3633         print *,'i',i,' sumprotype',sumpro_type(i)
3634       enddo
3635       return
3636       end subroutine mcmread
3637 !-----------------------------------------------------------------------------
3638       subroutine read_minim
3639
3640       use minim_data
3641 !      implicit real*8 (a-h,o-z)
3642 !      include 'DIMENSIONS'
3643 !      include 'COMMON.MINIM'
3644 !      include 'COMMON.IOUNITS'
3645 !      character(len=80) :: ucase
3646       character(len=320) :: minimcard
3647 !el local variables
3648 !     integer :: ntf,ik,iw_pdb
3649 !     real(kind=8) :: var,ene
3650
3651       call card_concat(minimcard,.true.)
3652       call readi(minimcard,'MAXMIN',maxmin,2000)
3653       call readi(minimcard,'MAXFUN',maxfun,5000)
3654       call readi(minimcard,'MINMIN',minmin,maxmin)
3655       call readi(minimcard,'MINFUN',minfun,maxmin)
3656       call reada(minimcard,'TOLF',tolf,1.0D-2)
3657       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3658       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3659       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3660       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3661       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3662                'Options in energy minimization:'
3663       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3664        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3665        'MinMin:',MinMin,' MinFun:',MinFun,&
3666        ' TolF:',TolF,' RTolF:',RTolF
3667       return
3668       end subroutine read_minim
3669 !-----------------------------------------------------------------------------
3670       subroutine openunits
3671
3672       use MD_data, only: usampl
3673       use csa_data
3674       use MPI_data
3675       use control_data, only:out1file
3676       use control, only: getenv_loc
3677 !      implicit real*8 (a-h,o-z)
3678 !      include 'DIMENSIONS'
3679 #ifdef MPI
3680       include 'mpif.h'
3681       character(len=16) :: form,nodename
3682       integer :: nodelen,ierror,npos
3683 #endif
3684 !      include 'COMMON.SETUP'
3685 !      include 'COMMON.IOUNITS'
3686 !      include 'COMMON.MD'
3687 !      include 'COMMON.CONTROL'
3688       integer :: lenpre,lenpot,lentmp   !,ilen
3689 !el      external ilen
3690       character(len=3) :: out1file_text !,ucase
3691       character(len=3) :: ll
3692 !el      external ucase
3693 !el local variables
3694 !     integer :: ntf,ik,iw_pdb
3695 !     real(kind=8) :: var,ene
3696 !
3697 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3698       call getenv_loc("PREFIX",prefix)
3699       pref_orig = prefix
3700       call getenv_loc("POT",pot)
3701       call getenv_loc("DIRTMP",tmpdir)
3702       call getenv_loc("CURDIR",curdir)
3703       call getenv_loc("OUT1FILE",out1file_text)
3704 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3705       out1file_text=ucase(out1file_text)
3706       if (out1file_text(1:1).eq."Y") then
3707         out1file=.true.
3708       else
3709         out1file=fg_rank.gt.0
3710       endif
3711       lenpre=ilen(prefix)
3712       lenpot=ilen(pot)
3713       lentmp=ilen(tmpdir)
3714       if (lentmp.gt.0) then
3715           write (*,'(80(1h!))')
3716           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
3717           write (*,'(80(1h!))')
3718           write (*,*)"All output files will be on node /tmp directory."
3719 #ifdef MPI
3720         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3721         if (me.eq.king) then
3722           write (*,*) "The master node is ",nodename
3723         else if (fg_rank.eq.0) then
3724           write (*,*) "I am the CG slave node ",nodename
3725         else
3726           write (*,*) "I am the FG slave node ",nodename
3727         endif
3728 #endif
3729         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3730         lenpre = lentmp+lenpre+1
3731       endif
3732       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3733 ! Get the names and open the input files
3734 #if defined(WINIFL) || defined(WINPGI)
3735       open(1,file=pref_orig(:ilen(pref_orig))// &
3736         '.inp',status='old',readonly,shared)
3737        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3738 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3739 ! Get parameter filenames and open the parameter files.
3740       call getenv_loc('BONDPAR',bondname)
3741       open (ibond,file=bondname,status='old',readonly,shared)
3742       call getenv_loc('THETPAR',thetname)
3743       open (ithep,file=thetname,status='old',readonly,shared)
3744       call getenv_loc('ROTPAR',rotname)
3745       open (irotam,file=rotname,status='old',readonly,shared)
3746       call getenv_loc('TORPAR',torname)
3747       open (itorp,file=torname,status='old',readonly,shared)
3748       call getenv_loc('TORDPAR',tordname)
3749       open (itordp,file=tordname,status='old',readonly,shared)
3750       call getenv_loc('FOURIER',fouriername)
3751       open (ifourier,file=fouriername,status='old',readonly,shared)
3752       call getenv_loc('ELEPAR',elename)
3753       open (ielep,file=elename,status='old',readonly,shared)
3754       call getenv_loc('SIDEPAR',sidename)
3755       open (isidep,file=sidename,status='old',readonly,shared)
3756 #elif (defined CRAY) || (defined AIX)
3757       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3758         action='read')
3759 !      print *,"Processor",myrank," opened file 1"
3760       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3761 !      print *,"Processor",myrank," opened file 9"
3762 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3763 ! Get parameter filenames and open the parameter files.
3764       call getenv_loc('BONDPAR',bondname)
3765       open (ibond,file=bondname,status='old',action='read')
3766 !      print *,"Processor",myrank," opened file IBOND"
3767       call getenv_loc('THETPAR',thetname)
3768       open (ithep,file=thetname,status='old',action='read')
3769 !      print *,"Processor",myrank," opened file ITHEP"
3770       call getenv_loc('ROTPAR',rotname)
3771       open (irotam,file=rotname,status='old',action='read')
3772 !      print *,"Processor",myrank," opened file IROTAM"
3773       call getenv_loc('TORPAR',torname)
3774       open (itorp,file=torname,status='old',action='read')
3775 !      print *,"Processor",myrank," opened file ITORP"
3776       call getenv_loc('TORDPAR',tordname)
3777       open (itordp,file=tordname,status='old',action='read')
3778 !      print *,"Processor",myrank," opened file ITORDP"
3779       call getenv_loc('SCCORPAR',sccorname)
3780       open (isccor,file=sccorname,status='old',action='read')
3781 !      print *,"Processor",myrank," opened file ISCCOR"
3782       call getenv_loc('FOURIER',fouriername)
3783       open (ifourier,file=fouriername,status='old',action='read')
3784 !      print *,"Processor",myrank," opened file IFOURIER"
3785       call getenv_loc('ELEPAR',elename)
3786       open (ielep,file=elename,status='old',action='read')
3787 !      print *,"Processor",myrank," opened file IELEP"
3788       call getenv_loc('SIDEPAR',sidename)
3789       open (isidep,file=sidename,status='old',action='read')
3790       call getenv_loc('LIPTRANPAR',liptranname)
3791       open (iliptranpar,file=liptranname,status='old',action='read')
3792       call getenv_loc('TUBEPAR',tubename)
3793       open (itube,file=tubename,status='old',action='read')
3794
3795 !      print *,"Processor",myrank," opened file ISIDEP"
3796 !      print *,"Processor",myrank," opened parameter files"
3797 #elif (defined G77)
3798       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3799       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3800 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3801 ! Get parameter filenames and open the parameter files.
3802       call getenv_loc('BONDPAR',bondname)
3803       open (ibond,file=bondname,status='old')
3804       call getenv_loc('THETPAR',thetname)
3805       open (ithep,file=thetname,status='old')
3806       call getenv_loc('ROTPAR',rotname)
3807       open (irotam,file=rotname,status='old')
3808       call getenv_loc('TORPAR',torname)
3809       open (itorp,file=torname,status='old')
3810       call getenv_loc('TORDPAR',tordname)
3811       open (itordp,file=tordname,status='old')
3812       call getenv_loc('SCCORPAR',sccorname)
3813       open (isccor,file=sccorname,status='old')
3814       call getenv_loc('FOURIER',fouriername)
3815       open (ifourier,file=fouriername,status='old')
3816       call getenv_loc('ELEPAR',elename)
3817       open (ielep,file=elename,status='old')
3818       call getenv_loc('SIDEPAR',sidename)
3819       open (isidep,file=sidename,status='old')
3820       call getenv_loc('LIPTRANPAR',liptranname)
3821       open (iliptranpar,file=liptranname,status='old')
3822       call getenv_loc('TUBEPAR',tubename)
3823       open (itube,file=tubename,status='old')
3824 #else
3825       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3826         readonly)
3827        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3828 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3829 ! Get parameter filenames and open the parameter files.
3830       call getenv_loc('BONDPAR',bondname)
3831       open (ibond,file=bondname,status='old',action='read')
3832       call getenv_loc('THETPAR',thetname)
3833       open (ithep,file=thetname,status='old',action='read')
3834       call getenv_loc('ROTPAR',rotname)
3835       open (irotam,file=rotname,status='old',action='read')
3836       call getenv_loc('TORPAR',torname)
3837       open (itorp,file=torname,status='old',action='read')
3838       call getenv_loc('TORDPAR',tordname)
3839       open (itordp,file=tordname,status='old',action='read')
3840       call getenv_loc('SCCORPAR',sccorname)
3841       open (isccor,file=sccorname,status='old',action='read')
3842 #ifndef CRYST_THETA
3843       call getenv_loc('THETPARPDB',thetname_pdb)
3844       print *,"thetname_pdb ",thetname_pdb
3845       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3846       print *,ithep_pdb," opened"
3847 #endif
3848       call getenv_loc('FOURIER',fouriername)
3849       open (ifourier,file=fouriername,status='old',readonly)
3850       call getenv_loc('ELEPAR',elename)
3851       open (ielep,file=elename,status='old',readonly)
3852       call getenv_loc('SIDEPAR',sidename)
3853       open (isidep,file=sidename,status='old',readonly)
3854       call getenv_loc('LIPTRANPAR',liptranname)
3855       open (iliptranpar,file=liptranname,status='old',action='read')
3856       call getenv_loc('TUBEPAR',tubename)
3857       open (itube,file=tubename,status='old',action='read')
3858
3859 #ifndef CRYST_SC
3860       call getenv_loc('ROTPARPDB',rotname_pdb)
3861       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3862 #endif
3863 #endif
3864 #ifndef OLDSCP
3865 !
3866 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3867 ! Use -DOLDSCP to use hard-coded constants instead.
3868 !
3869       call getenv_loc('SCPPAR',scpname)
3870 #if defined(WINIFL) || defined(WINPGI)
3871       open (iscpp,file=scpname,status='old',readonly,shared)
3872 #elif (defined CRAY)  || (defined AIX)
3873       open (iscpp,file=scpname,status='old',action='read')
3874 #elif (defined G77)
3875       open (iscpp,file=scpname,status='old')
3876 #else
3877       open (iscpp,file=scpname,status='old',action='read')
3878 #endif
3879 #endif
3880       call getenv_loc('PATTERN',patname)
3881 #if defined(WINIFL) || defined(WINPGI)
3882       open (icbase,file=patname,status='old',readonly,shared)
3883 #elif (defined CRAY)  || (defined AIX)
3884       open (icbase,file=patname,status='old',action='read')
3885 #elif (defined G77)
3886       open (icbase,file=patname,status='old')
3887 #else
3888       open (icbase,file=patname,status='old',action='read')
3889 #endif
3890 #ifdef MPI
3891 ! Open output file only for CG processes
3892 !      print *,"Processor",myrank," fg_rank",fg_rank
3893       if (fg_rank.eq.0) then
3894
3895       if (nodes.eq.1) then
3896         npos=3
3897       else
3898         npos = dlog10(dfloat(nodes-1))+1
3899       endif
3900       if (npos.lt.3) npos=3
3901       write (liczba,'(i1)') npos
3902       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3903         //')'
3904       write (liczba,form) me
3905       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3906         liczba(:ilen(liczba))
3907       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3908         //'.int'
3909       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3910         //'.pdb'
3911       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3912         liczba(:ilen(liczba))//'.mol2'
3913       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3914         liczba(:ilen(liczba))//'.stat'
3915       if (lentmp.gt.0) &
3916         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3917             //liczba(:ilen(liczba))//'.stat')
3918       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3919         //'.rst'
3920       if(usampl) then
3921           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3922        liczba(:ilen(liczba))//'.const'
3923       endif
3924
3925       endif
3926 #else
3927       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3928       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3929       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3930       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3931       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3932       if (lentmp.gt.0) &
3933         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3934           '.stat')
3935       rest2name=prefix(:ilen(prefix))//'.rst'
3936       if(usampl) then
3937          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3938       endif
3939 #endif
3940 #if defined(AIX) || defined(PGI)
3941       if (me.eq.king .or. .not. out1file) &
3942          open(iout,file=outname,status='unknown')
3943 #ifdef DEBUG
3944       if (fg_rank.gt.0) then
3945         write (liczba,'(i3.3)') myrank/nfgtasks
3946         write (ll,'(bz,i3.3)') fg_rank
3947         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3948          status='unknown')
3949       endif
3950 #endif
3951       if(me.eq.king) then
3952        open(igeom,file=intname,status='unknown',position='append')
3953        open(ipdb,file=pdbname,status='unknown')
3954        open(imol2,file=mol2name,status='unknown')
3955        open(istat,file=statname,status='unknown',position='append')
3956       else
3957 !1out       open(iout,file=outname,status='unknown')
3958       endif
3959 #else
3960       if (me.eq.king .or. .not.out1file) &
3961           open(iout,file=outname,status='unknown')
3962 #ifdef DEBUG
3963       if (fg_rank.gt.0) then
3964         write (liczba,'(i3.3)') myrank/nfgtasks
3965         write (ll,'(bz,i3.3)') fg_rank
3966         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3967          status='unknown')
3968       endif
3969 #endif
3970       if(me.eq.king) then
3971        open(igeom,file=intname,status='unknown',access='append')
3972        open(ipdb,file=pdbname,status='unknown')
3973        open(imol2,file=mol2name,status='unknown')
3974        open(istat,file=statname,status='unknown',access='append')
3975       else
3976 !1out       open(iout,file=outname,status='unknown')
3977       endif
3978 #endif
3979       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3980       csa_seed=prefix(:lenpre)//'.CSA.seed'
3981       csa_history=prefix(:lenpre)//'.CSA.history'
3982       csa_bank=prefix(:lenpre)//'.CSA.bank'
3983       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3984       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3985       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3986 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3987       csa_int=prefix(:lenpre)//'.int'
3988       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3989       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3990       csa_in=prefix(:lenpre)//'.CSA.in'
3991 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3992 ! Write file names
3993       if (me.eq.king)then
3994       write (iout,'(80(1h-))')
3995       write (iout,'(30x,a)') "FILE ASSIGNMENT"
3996       write (iout,'(80(1h-))')
3997       write (iout,*) "Input file                      : ",&
3998         pref_orig(:ilen(pref_orig))//'.inp'
3999       write (iout,*) "Output file                     : ",&
4000         outname(:ilen(outname))
4001       write (iout,*)
4002       write (iout,*) "Sidechain potential file        : ",&
4003         sidename(:ilen(sidename))
4004 #ifndef OLDSCP
4005       write (iout,*) "SCp potential file              : ",&
4006         scpname(:ilen(scpname))
4007 #endif
4008       write (iout,*) "Electrostatic potential file    : ",&
4009         elename(:ilen(elename))
4010       write (iout,*) "Cumulant coefficient file       : ",&
4011         fouriername(:ilen(fouriername))
4012       write (iout,*) "Torsional parameter file        : ",&
4013         torname(:ilen(torname))
4014       write (iout,*) "Double torsional parameter file : ",&
4015         tordname(:ilen(tordname))
4016       write (iout,*) "SCCOR parameter file : ",&
4017         sccorname(:ilen(sccorname))
4018       write (iout,*) "Bond & inertia constant file    : ",&
4019         bondname(:ilen(bondname))
4020       write (iout,*) "Bending parameter file          : ",&
4021         thetname(:ilen(thetname))
4022       write (iout,*) "Rotamer parameter file          : ",&
4023         rotname(:ilen(rotname))
4024 !el----
4025 #ifndef CRYST_THETA
4026       write (iout,*) "Thetpdb parameter file          : ",&
4027         thetname_pdb(:ilen(thetname_pdb))
4028 #endif
4029 !el
4030       write (iout,*) "Threading database              : ",&
4031         patname(:ilen(patname))
4032       if (lentmp.ne.0) &
4033       write (iout,*)" DIRTMP                          : ",&
4034         tmpdir(:lentmp)
4035       write (iout,'(80(1h-))')
4036       endif
4037       return
4038       end subroutine openunits
4039 !-----------------------------------------------------------------------------
4040       subroutine readrst
4041
4042       use geometry_data, only: nres,dc
4043       use MD_data
4044 !      implicit real*8 (a-h,o-z)
4045 !      include 'DIMENSIONS'
4046 !      include 'COMMON.CHAIN'
4047 !      include 'COMMON.IOUNITS'
4048 !      include 'COMMON.MD'
4049 !el local variables
4050       integer ::i,j
4051 !     real(kind=8) :: var,ene
4052
4053       open(irest2,file=rest2name,status='unknown')
4054       read(irest2,*) totT,EK,potE,totE,t_bath
4055 !      do i=1,2*nres
4056 ! AL 4/17/17: Now reading d_t(0,:) too
4057       do i=0,2*nres
4058          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4059       enddo
4060 !      do i=1,2*nres
4061 ! AL 4/17/17: Now reading d_c(0,:) too
4062       do i=0,2*nres
4063          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4064       enddo
4065       if(usampl) then
4066              read (irest2,*) iset
4067       endif
4068       close(irest2)
4069       return
4070       end subroutine readrst
4071 !-----------------------------------------------------------------------------
4072       subroutine read_fragments
4073
4074       use energy_data
4075 !      use geometry
4076       use control_data, only:out1file
4077       use MD_data
4078       use MPI_data
4079 !      implicit real*8 (a-h,o-z)
4080 !      include 'DIMENSIONS'
4081 #ifdef MPI
4082       include 'mpif.h'
4083 #endif
4084 !      include 'COMMON.SETUP'
4085 !      include 'COMMON.CHAIN'
4086 !      include 'COMMON.IOUNITS'
4087 !      include 'COMMON.MD'
4088 !      include 'COMMON.CONTROL'
4089 !el local variables
4090       integer :: i
4091 !     real(kind=8) :: var,ene
4092
4093       read(inp,*) nset,nfrag,npair,nfrag_back
4094
4095 !el from module energy
4096 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
4097       if(.not.allocated(wfrag_back)) then
4098         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4099         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4100
4101         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4102         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
4103
4104         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4105         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
4106       endif
4107
4108       if(me.eq.king.or..not.out1file) &
4109        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4110         " nfrag_back",nfrag_back
4111       do iset=1,nset
4112          read(inp,*) mset(iset)
4113        do i=1,nfrag
4114          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4115            qinfrag(i,iset)
4116          if(me.eq.king.or..not.out1file) &
4117           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4118            ifrag(2,i,iset), qinfrag(i,iset)
4119        enddo
4120        do i=1,npair
4121         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4122           qinpair(i,iset)
4123         if(me.eq.king.or..not.out1file) &
4124          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4125           ipair(2,i,iset), qinpair(i,iset)
4126        enddo
4127        do i=1,nfrag_back
4128         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4129            wfrag_back(3,i,iset),&
4130            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4131         if(me.eq.king.or..not.out1file) &
4132          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4133          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4134        enddo
4135       enddo
4136       return
4137       end subroutine read_fragments
4138 !-----------------------------------------------------------------------------
4139 ! shift.F       io_csa
4140 !-----------------------------------------------------------------------------
4141       subroutine csa_read
4142
4143       use csa_data
4144 !      implicit real*8 (a-h,o-z)
4145 !      include 'DIMENSIONS'
4146 !      include 'COMMON.CSA'
4147 !      include 'COMMON.BANK'
4148 !      include 'COMMON.IOUNITS'
4149 !el local variables
4150 !     integer :: ntf,ik,iw_pdb
4151 !     real(kind=8) :: var,ene
4152
4153       open(icsa_in,file=csa_in,status="old",err=100)
4154        read(icsa_in,*) nconf
4155        read(icsa_in,*) jstart,jend
4156        read(icsa_in,*) nstmax
4157        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4158        read(icsa_in,*) nran0,nran1,irr
4159        read(icsa_in,*) nseed
4160        read(icsa_in,*) ntotal,cut1,cut2
4161        read(icsa_in,*) estop
4162        read(icsa_in,*) icmax,irestart
4163        read(icsa_in,*) ntbankm,dele,difcut
4164        read(icsa_in,*) iref,rmscut,pnccut
4165        read(icsa_in,*) ndiff
4166       close(icsa_in)
4167
4168       return
4169
4170  100  continue
4171       return
4172       end subroutine csa_read
4173 !-----------------------------------------------------------------------------
4174       subroutine initial_write
4175
4176       use csa_data
4177 !      implicit real*8 (a-h,o-z)
4178 !      include 'DIMENSIONS'
4179 !      include 'COMMON.CSA'
4180 !      include 'COMMON.BANK'
4181 !      include 'COMMON.IOUNITS'
4182 !el local variables
4183 !     integer :: ntf,ik,iw_pdb
4184 !     real(kind=8) :: var,ene
4185
4186       open(icsa_seed,file=csa_seed,status="unknown")
4187        write(icsa_seed,*) "seed"
4188       close(31)
4189 #if defined(AIX) || defined(PGI)
4190        open(icsa_history,file=csa_history,status="unknown",&
4191         position="append")
4192 #else
4193        open(icsa_history,file=csa_history,status="unknown",&
4194         access="append")
4195 #endif
4196        write(icsa_history,*) nconf
4197        write(icsa_history,*) jstart,jend
4198        write(icsa_history,*) nstmax
4199        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4200        write(icsa_history,*) nran0,nran1,irr
4201        write(icsa_history,*) nseed
4202        write(icsa_history,*) ntotal,cut1,cut2
4203        write(icsa_history,*) estop
4204        write(icsa_history,*) icmax,irestart
4205        write(icsa_history,*) ntbankm,dele,difcut
4206        write(icsa_history,*) iref,rmscut,pnccut
4207        write(icsa_history,*) ndiff
4208
4209        write(icsa_history,*)
4210       close(icsa_history)
4211
4212       open(icsa_bank1,file=csa_bank1,status="unknown")
4213        write(icsa_bank1,*) 0
4214       close(icsa_bank1)
4215
4216       return
4217       end subroutine initial_write
4218 !-----------------------------------------------------------------------------
4219       subroutine restart_write
4220
4221       use csa_data
4222 !      implicit real*8 (a-h,o-z)
4223 !      include 'DIMENSIONS'
4224 !      include 'COMMON.IOUNITS'
4225 !      include 'COMMON.CSA'
4226 !      include 'COMMON.BANK'
4227 !el local variables
4228 !     integer :: ntf,ik,iw_pdb
4229 !     real(kind=8) :: var,ene
4230
4231 #if defined(AIX) || defined(PGI)
4232        open(icsa_history,file=csa_history,position="append")
4233 #else
4234        open(icsa_history,file=csa_history,access="append")
4235 #endif
4236        write(icsa_history,*)
4237        write(icsa_history,*) "This is restart"
4238        write(icsa_history,*)
4239        write(icsa_history,*) nconf
4240        write(icsa_history,*) jstart,jend
4241        write(icsa_history,*) nstmax
4242        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4243        write(icsa_history,*) nran0,nran1,irr
4244        write(icsa_history,*) nseed
4245        write(icsa_history,*) ntotal,cut1,cut2
4246        write(icsa_history,*) estop
4247        write(icsa_history,*) icmax,irestart
4248        write(icsa_history,*) ntbankm,dele,difcut
4249        write(icsa_history,*) iref,rmscut,pnccut
4250        write(icsa_history,*) ndiff
4251        write(icsa_history,*)
4252        write(icsa_history,*) "irestart is: ", irestart
4253
4254        write(icsa_history,*)
4255       close(icsa_history)
4256
4257       return
4258       end subroutine restart_write
4259 !-----------------------------------------------------------------------------
4260 ! test.F
4261 !-----------------------------------------------------------------------------
4262       subroutine write_pdb(npdb,titelloc,ee)
4263
4264 !      implicit real*8 (a-h,o-z)
4265 !      include 'DIMENSIONS'
4266 !      include 'COMMON.IOUNITS'
4267       character(len=50) :: titelloc1
4268       character*(*) :: titelloc
4269       character(len=3) :: zahl
4270       character(len=5) :: liczba5
4271       real(kind=8) :: ee
4272       integer :: npdb   !,ilen
4273 !el      external ilen
4274 !el local variables
4275       integer :: lenpre
4276 !     real(kind=8) :: var,ene
4277
4278       titelloc1=titelloc
4279       lenpre=ilen(prefix)
4280       if (npdb.lt.1000) then
4281        call numstr(npdb,zahl)
4282        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4283       else
4284         if (npdb.lt.10000) then
4285          write(liczba5,'(i1,i4)') 0,npdb
4286         else
4287          write(liczba5,'(i5)') npdb
4288         endif
4289         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4290       endif
4291       call pdbout(ee,titelloc1,ipdb)
4292       close(ipdb)
4293       return
4294       end subroutine write_pdb
4295 !-----------------------------------------------------------------------------
4296 ! thread.F
4297 !-----------------------------------------------------------------------------
4298       subroutine write_thread_summary
4299 ! Thread the sequence through a database of known structures
4300       use control_data, only: refstr
4301 !      use geometry
4302       use energy_data, only: n_ene_comp
4303       use compare_data
4304 !      implicit real*8 (a-h,o-z)
4305 !      include 'DIMENSIONS'
4306 #ifdef MPI
4307       use MPI_data      !include 'COMMON.INFO'
4308       include 'mpif.h'
4309 #endif
4310 !      include 'COMMON.CONTROL'
4311 !      include 'COMMON.CHAIN'
4312 !      include 'COMMON.DBASE'
4313 !      include 'COMMON.INTERACT'
4314 !      include 'COMMON.VAR'
4315 !      include 'COMMON.THREAD'
4316 !      include 'COMMON.FFIELD'
4317 !      include 'COMMON.SBRIDGE'
4318 !      include 'COMMON.HEADER'
4319 !      include 'COMMON.NAMES'
4320 !      include 'COMMON.IOUNITS'
4321 !      include 'COMMON.TIME1'
4322
4323       integer,dimension(maxthread) :: ip
4324       real(kind=8),dimension(0:n_ene) :: energia
4325 !el local variables
4326       integer :: i,j,ii,jj,ipj,ik,kk,ist
4327       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4328
4329       write (iout,'(30x,a/)') &
4330        '  *********** Summary threading statistics ************'
4331       write (iout,'(a)') 'Initial energies:'
4332       write (iout,'(a4,2x,a12,14a14,3a8)') &
4333        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4334        'RMSnat','NatCONT','NNCONT','RMS'
4335 ! Energy sort patterns
4336       do i=1,nthread
4337         ip(i)=i
4338       enddo
4339       do i=1,nthread-1
4340         enet=ener(n_ene-1,ip(i))
4341         jj=i
4342         do j=i+1,nthread
4343           if (ener(n_ene-1,ip(j)).lt.enet) then
4344             jj=j
4345             enet=ener(n_ene-1,ip(j))
4346           endif
4347         enddo
4348         if (jj.ne.i) then
4349           ipj=ip(jj)
4350           ip(jj)=ip(i)
4351           ip(i)=ipj
4352         endif
4353       enddo
4354       do ik=1,nthread
4355         i=ip(ik)
4356         ii=ipatt(1,i)
4357         ist=nres_base(2,ii)+ipatt(2,i)
4358         do kk=1,n_ene_comp
4359           energia(i)=ener0(kk,i)
4360         enddo
4361         etot=ener0(n_ene_comp+1,i)
4362         rmsnat=ener0(n_ene_comp+2,i)
4363         rms=ener0(n_ene_comp+3,i)
4364         frac=ener0(n_ene_comp+4,i)
4365         frac_nn=ener0(n_ene_comp+5,i)
4366
4367         if (refstr) then
4368         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4369         i,str_nam(ii),ist+1,&
4370         (energia(print_order(kk)),kk=1,nprint_ene),&
4371         etot,rmsnat,frac,frac_nn,rms
4372         else
4373         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4374         i,str_nam(ii),ist+1,&
4375         (energia(print_order(kk)),kk=1,nprint_ene),etot
4376         endif
4377       enddo
4378       write (iout,'(//a)') 'Final energies:'
4379       write (iout,'(a4,2x,a12,17a14,3a8)') &
4380        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4381        'RMSnat','NatCONT','NNCONT','RMS'
4382       do ik=1,nthread
4383         i=ip(ik)
4384         ii=ipatt(1,i)
4385         ist=nres_base(2,ii)+ipatt(2,i)
4386         do kk=1,n_ene_comp
4387           energia(kk)=ener(kk,ik)
4388         enddo
4389         etot=ener(n_ene_comp+1,i)
4390         rmsnat=ener(n_ene_comp+2,i)
4391         rms=ener(n_ene_comp+3,i)
4392         frac=ener(n_ene_comp+4,i)
4393         frac_nn=ener(n_ene_comp+5,i)
4394         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4395         i,str_nam(ii),ist+1,&
4396         (energia(print_order(kk)),kk=1,nprint_ene),&
4397         etot,rmsnat,frac,frac_nn,rms
4398       enddo
4399       write (iout,'(/a/)') 'IEXAM array:'
4400       write (iout,'(i5)') nexcl
4401       do i=1,nexcl
4402         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4403       enddo
4404       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4405        'Max. time for threading step ',max_time_for_thread,&
4406        'Average time for threading step: ',ave_time_for_thread
4407       return
4408       end subroutine write_thread_summary
4409 !-----------------------------------------------------------------------------
4410       subroutine write_stat_thread(ithread,ipattern,ist)
4411
4412       use energy_data, only: n_ene_comp
4413       use compare_data
4414 !      implicit real*8 (a-h,o-z)
4415 !      include "DIMENSIONS"
4416 !      include "COMMON.CONTROL"
4417 !      include "COMMON.IOUNITS"
4418 !      include "COMMON.THREAD"
4419 !      include "COMMON.FFIELD"
4420 !      include "COMMON.DBASE"
4421 !      include "COMMON.NAMES"
4422       real(kind=8),dimension(0:n_ene) :: energia
4423 !el local variables
4424       integer :: ithread,ipattern,ist,i
4425       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4426
4427 #if defined(AIX) || defined(PGI)
4428       open(istat,file=statname,position='append')
4429 #else
4430       open(istat,file=statname,access='append')
4431 #endif
4432       do i=1,n_ene_comp
4433         energia(i)=ener(i,ithread)
4434       enddo
4435       etot=ener(n_ene_comp+1,ithread)
4436       rmsnat=ener(n_ene_comp+2,ithread)
4437       rms=ener(n_ene_comp+3,ithread)
4438       frac=ener(n_ene_comp+4,ithread)
4439       frac_nn=ener(n_ene_comp+5,ithread)
4440       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4441         ithread,str_nam(ipattern),ist+1,&
4442         (energia(print_order(i)),i=1,nprint_ene),&
4443         etot,rmsnat,frac,frac_nn,rms
4444       close (istat)
4445       return
4446       end subroutine write_stat_thread
4447 !-----------------------------------------------------------------------------
4448 #endif
4449 !-----------------------------------------------------------------------------
4450       end module io_config