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