source/unres/io_config.f90

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