source/unres/io_config.f90

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