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 !defined(WHAM_RUN) && !defined(CLUSTER)
  34 !-----------------------------------------------------------------------------
  35 ! bank.F    io_csa
  36 !-----------------------------------------------------------------------------
  37       subroutine write_rbank(jlee,adif,nft)
  38
  39       use csa_data
  40       use geometry_data, only: nres,rad2deg
  41 !      implicit real*8 (a-h,o-z)
  42 !      include 'DIMENSIONS'
  43 !      include 'COMMON.IOUNITS'
  44 !      include 'COMMON.CSA'
  45 !      include 'COMMON.BANK'
  46 !      include 'COMMON.CHAIN'
  47 !      include 'COMMON.GEO'
  48 !el local variables
  49       integer :: nft,i,k,j,l,jlee
  50       real(kind=8) :: adif
  51
  52       open(icsa_rbank,file=csa_rbank,status="unknown")
  53       write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
  54       do k=1,nbank
  55        write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
  56        do j=1,numch
  57         do l=2,nres-1
  58          write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
  59         enddo
  60        enddo
  61       enddo
  62       close(icsa_rbank)
  63
  64   850 format (10f8.3)
  65   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
  66               i8,i10,i2,f15.5)
  67   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
  68                   ' %NC ',0pf5.2)
  69
  70       return
  71       end subroutine write_rbank
  72 !-----------------------------------------------------------------------------
  73       subroutine read_rbank(jlee,adif)
  74
  75       use csa_data
  76       use geometry_data, only: nres,deg2rad
  77       use MPI_data
  78 !      implicit real*8 (a-h,o-z)
  79 !      include 'DIMENSIONS'
  80       include 'mpif.h'
  81 !      include 'COMMON.IOUNITS'
  82 !      include 'COMMON.CSA'
  83 !      include 'COMMON.BANK'
  84 !      include 'COMMON.CHAIN'
  85 !      include 'COMMON.GEO'
  86 !      include 'COMMON.SETUP'
  87       character(len=80) :: karta
  88 !el local variables
  89       integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
  90                  ierror,ierrcode,jlee,jleer
  91       real(kind=8) :: adif
  92
  93       open(icsa_rbank,file=csa_rbank,status="old")
  94       read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
  95       print *,jleer,nbankr,nstepr,nftr,icycler,adif
  96 !       print *, 'adif from read_rbank ',adif
  97 #ifdef MPI
  98       if(nbankr.ne.nbank) then
  99         write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
 100         ' NBANK',nbank
 101         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 102       endif
 103       if(jleer.ne.jlee) then
 104         write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
 105         ' JLEE',jlee
 106         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 107       endif
 108 #endif
 109
 110       kk=0
 111       do k=1,nbankr
 112         read (icsa_rbank,'(a80)') karta
 113         write(iout,*) "READ_RBANK: kk=",kk
 114         write(iout,*) karta
 115 !        if (index(karta,"*").gt.0) then
 116 !          write (iout,*) "***** Stars in bankr ***** k=",k,
 117 !     &      " skipped"
 118 !          do j=1,numch
 119 !            do l=2,nres-1
 120 !              read (30,850) (rdummy,i=1,4)
 121 !            enddo
 122 !          enddo
 123 !        else
 124           kk=kk+1
 125           call reada(karta,"total E",rene(kk),1.0d20)
 126           call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
 127           call reada(karta,"%NC",rpncn(kk),0.0d0)
 128           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
 129             "%NC",bpncn(kk),ibank(kk)
 130 !          read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
 131           do j=1,numch
 132             do l=2,nres-1
 133               read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
 134 !              write (iout,850) (rvar(i,l,j,kk),i=1,4)
 135               do i=1,4
 136                 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
 137               enddo
 138             enddo
 139           enddo
 140 !        endif
 141       enddo
 142 !d      write (*,*) "read_rbank ******************* kk",kk,
 143 !d     &  "nbankr",nbankr
 144       if (kk.lt.nbankr) nbankr=kk
 145 !d      do kk=1,nbankr
 146 !d          print *,"kk=",kk
 147 !d          do j=1,numch
 148 !d            do l=2,nres-1
 149 !d              write (*,850) (rvar(i,l,j,kk),i=1,4)
 150 !d            enddo
 151 !d          enddo
 152 !d      enddo
 153       close(icsa_rbank)
 154
 155   850 format (10f8.3)
 156   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 157   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
 158
 159       return
 160       end subroutine read_rbank
 161 !-----------------------------------------------------------------------------
 162       subroutine write_bank(jlee,nft)
 163
 164       use csa_data
 165       use control_data, only: vdisulf
 166       use geometry_data, only: nres,rad2deg
 167 !      implicit real*8 (a-h,o-z)
 168 !      include 'DIMENSIONS'
 169 !      include 'COMMON.IOUNITS'
 170 !      include 'COMMON.CSA'
 171 !      include 'COMMON.BANK'
 172 !      include 'COMMON.CHAIN'
 173 !      include 'COMMON.GEO'
 174 !      include 'COMMON.SBRIDGE'
 175 !     include 'COMMON.CONTROL'
 176       character(len=7) :: chtmp
 177       character(len=40) :: chfrm
 178 !el      external ilen
 179 !el local variables
 180       integer :: nft,k,l,i,j,jlee
 181
 182       open(icsa_bank,file=csa_bank,status="unknown")
 183       write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
 184       write (icsa_bank,902) nglob_csa, eglob_csa
 185       open (igeom,file=intname,status='UNKNOWN')
 186       do k=1,nbank
 187        write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 188        if (vdisulf) write (icsa_bank,'(101i4)') &
 189           bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
 190        do j=1,numch
 191         do l=2,nres-1
 192          write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 193         enddo
 194        enddo
 195        if (bvar_nss(k).le.9) then
 196          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 197            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 198        else
 199          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 200            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 201          write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
 202                                       bvar_ss(2,i,k),i=10,bvar_nss(k))
 203        endif
 204        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 205        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 206        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 207        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 208       enddo
 209       close(icsa_bank)
 210       close(igeom)
 211
 212       if (nstep/200.gt.ilastnstep) then
 213
 214        ilastnstep=(ilastnstep+1)*1.5
 215        write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
 216        write(chtmp,chfrm) nstep
 217        open(icsa_int,file=prefix(:ilen(prefix)) &
 218                //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
 219        do k=1,nbank
 220         if (bvar_nss(k).le.9) then
 221          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 222         bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 223         else
 224          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 225            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 226          write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
 227                                 bvar_ss(2,i,k),i=10,bvar_nss(k))
 228         endif
 229         write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 230         write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 231         write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 232         write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 233        enddo
 234        close(icsa_int)
 235       endif
 236
 237
 238   200 format (8f10.4)
 239   850 format (10f8.3)
 240   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
 241               i8,i10,i2,f15.5)
 242   902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
 243   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 244               ' %NC ',0pf5.2,i5)
 245
 246       return
 247       end subroutine write_bank
 248 !-----------------------------------------------------------------------------
 249       subroutine write_bank_reminimized(jlee,nft)
 250
 251       use csa_data
 252       use geometry_data, only: nres,rad2deg
 253       use energy_data
 254 !      implicit real*8 (a-h,o-z)
 255 !      include 'DIMENSIONS'
 256 !      include 'COMMON.IOUNITS'
 257 !      include 'COMMON.CSA'
 258 !      include 'COMMON.BANK'
 259 !      include 'COMMON.CHAIN'
 260 !      include 'COMMON.GEO'
 261 !      include 'COMMON.SBRIDGE'
 262 !el local variables
 263       integer :: nft,i,l,j,k,jlee
 264
 265       open(icsa_bank_reminimized,file=csa_bank_reminimized,&
 266         status="unknown")
 267       write (icsa_bank_reminimized,900) &
 268         jlee,nbank,nstep,nft,icycle,cutdif
 269       open (igeom,file=intname,status='UNKNOWN')
 270       do k=1,nbank
 271        write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
 272         bpncn(k),ibank(k)
 273        do j=1,numch
 274         do l=2,nres-1
 275          write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 276         enddo
 277        enddo
 278        if (nss.le.9) then
 279          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 280            nss,(ihpb(i),jhpb(i),i=1,nss)
 281        else
 282          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 283            nss,(ihpb(i),jhpb(i),i=1,9)
 284          write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
 285        endif
 286        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 287        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 288        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 289        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 290       enddo
 291       close(icsa_bank_reminimized)
 292       close(igeom)
 293
 294   200 format (8f10.4)
 295   850 format (10f8.3)
 296   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
 297               i8,i10,i2,f15.5)
 298   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 299                ' %NC ',0pf5.2,i5)
 300
 301       return
 302       end subroutine write_bank_reminimized
 303 !-----------------------------------------------------------------------------
 304       subroutine read_bank(jlee,nft,cutdifr)
 305
 306       use csa_data
 307       use control_data, only: vdisulf
 308       use geometry_data, only: nres,deg2rad
 309       use energy_data
 310 !      implicit real*8 (a-h,o-z)
 311 !      include 'DIMENSIONS'
 312 !      include 'COMMON.IOUNITS'
 313 !      include 'COMMON.CSA'
 314 !      include 'COMMON.BANK'
 315 !      include 'COMMON.CHAIN'
 316 !      include 'COMMON.GEO'
 317 !      include 'COMMON.CONTROL'
 318 !      include 'COMMON.SBRIDGE'
 319       character(len=80) :: karta
 320 !      integer ilen
 321 !el      external ilen
 322 !el local variables
 323       integer :: nft,kk,k,l,i,j,jlee
 324       real(kind=8) :: cutdifr
 325
 326       open(icsa_bank,file=csa_bank,status="old")
 327        read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
 328        read (icsa_bank,902) nglob_csa, eglob_csa
 329 !      if(jleer.ne.jlee) then
 330 !        write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
 331 !    &   ' JLEE',jlee
 332 !        call mpi_abort(mpi_comm_world,ierror,ierrcode)
 333 !      endif
 334
 335       kk=0
 336       do k=1,nbank
 337         read (icsa_bank,'(a80)') karta
 338         write(iout,*) "READ_BANK: kk=",kk
 339         write(iout,*) karta
 340 !        if (index(karta,"*").gt.0) then
 341 !          write (iout,*) "***** Stars in bank ***** k=",k,
 342 !     &      " skipped"
 343 !          do j=1,numch
 344 !            do l=2,nres-1
 345 !              read (33,850) (rdummy,i=1,4)
 346 !            enddo
 347 !          enddo
 348 !        else
 349           kk=kk+1
 350           call reada(karta,"total E",bene(kk),1.0d20)
 351           call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
 352           call reada(karta,"%NC",bpncn(kk),0.0d0)
 353           read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
 354           goto 112
 355   111     ibank(kk)=0
 356   112     continue
 357           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
 358             "%NC",bpncn(kk),ibank(kk)
 359 !          read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 360           if (vdisulf) then
 361             read (icsa_bank,'(101i4)') &
 362               bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 363             bvar_ns(kk)=ns-2*bvar_nss(kk)
 364             write(iout,*) 'read SSBOND',bvar_nss(kk),&
 365                           ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 366 !d          write(iout,*) 'read CYS #free ', bvar_ns(kk)
 367             l=0
 368             do i=1,ns
 369              j=1
 370              do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
 371                        iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
 372                        j.le.bvar_nss(kk))
 373                 j=j+1
 374              enddo
 375              if (j.gt.bvar_nss(kk)) then
 376                l=l+1
 377                bvar_s(l,kk)=iss(i)
 378              endif
 379             enddo
 380 !d            write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
 381           endif
 382           do j=1,numch
 383             do l=2,nres-1
 384               read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
 385 !              write (iout,850) (bvar(i,l,j,kk),i=1,4)
 386               do i=1,4
 387                 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
 388               enddo ! l
 389             enddo ! l
 390           enddo ! j
 391 !        endif
 392       enddo ! k
 393
 394       if (kk.lt.nbank) nbank=kk
 395 !d      write (*,*) "read_bank ******************* kk",kk,
 396 !d     &  "nbank",nbank
 397 !d      do kk=1,nbank
 398 !d          print *,"kk=",kk
 399 !d          do j=1,numch
 400 !d            do l=2,nres-1
 401 !d              write (*,850) (bvar(i,l,j,kk),i=1,4)
 402 !d            enddo
 403 !d          enddo
 404 !d      enddo
 405
 406 !       do k=1,nbank
 407 !        read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 408 !        do j=1,numch
 409 !         do l=2,nres-1
 410 !          read (33,850) (bvar(i,l,j,k),i=1,4)
 411 !          do i=1,4
 412 !           bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
 413 !          enddo
 414 !         enddo
 415 !        enddo
 416 !       enddo
 417       close(icsa_bank)
 418
 419   850 format (10f8.3)
 420   952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
 421   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 422   902 format (1x,11x,i4,12x,1pe14.5)
 423   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
 424
 425       return
 426       end subroutine read_bank
 427 !-----------------------------------------------------------------------------
 428       subroutine write_bank1(jlee)
 429
 430       use csa_data
 431       use geometry_data, only: nres,rad2deg
 432 !      implicit real*8 (a-h,o-z)
 433 !      include 'DIMENSIONS'
 434 !      include 'COMMON.IOUNITS'
 435 !      include 'COMMON.CSA'
 436 !      include 'COMMON.BANK'
 437 !      include 'COMMON.CHAIN'
 438 !      include 'COMMON.GEO'
 439 !el local variables
 440       integer :: k,i,l,j,jlee
 441
 442 #if defined(AIX) || defined(PGI)
 443       open(icsa_bank1,file=csa_bank1,position="append")
 444 #else
 445       open(icsa_bank1,file=csa_bank1,access="append")
 446 #endif
 447       write (icsa_bank1,900) jlee,nbank,nstep,cutdif
 448       do k=1,nbank
 449        write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 450        do j=1,numch
 451         do l=2,nres-1
 452          write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 453         enddo
 454        enddo
 455       enddo
 456       close(icsa_bank1)
 457   850 format (10f8.3)
 458   900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
 459   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 460               ' %NC ',0pf5.2,i5)
 461
 462       return
 463       end subroutine write_bank1
 464 !-----------------------------------------------------------------------------
 465 ! cartprint.f
 466 !-----------------------------------------------------------------------------
 467 !      subroutine cartprint
 468
 469 !      use geometry_data, only: c
 470 !      use energy_data, only: itype
 471 !      implicit real*8 (a-h,o-z)
 472 !      include 'DIMENSIONS'
 473 !      include 'COMMON.CHAIN'
 474 !      include 'COMMON.INTERACT'
 475 !      include 'COMMON.NAMES'
 476 !      include 'COMMON.IOUNITS'
 477 !      integer :: i
 478
 479 !     write (iout,100)
 480 !      do i=1,nres
 481 !        write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
 482 !          c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
 483 !      enddo
 484 !  100 format (//'              alpha-carbon coordinates       ',&
 485 !                '     centroid coordinates'/ &
 486 !                '       ', 6X,'X',11X,'Y',11X,'Z',&
 487 !                                10X,'X',11X,'Y',11X,'Z')
 488 !  110 format (a,'(',i3,')',6f12.5)
 489 !      return
 490 !      end subroutine cartprint
 491 !-----------------------------------------------------------------------------
 492 ! dihed_cons.F
 493 !-----------------------------------------------------------------------------
 494       subroutine secstrp2dihc
 495
 496       use geometry_data
 497       use energy_data
 498 !      implicit real*8 (a-h,o-z)
 499 !      include 'DIMENSIONS'
 500 !      include 'COMMON.GEO'
 501 !      include 'COMMON.BOUNDS'
 502 !      include 'COMMON.CHAIN'
 503 !      include 'COMMON.TORCNSTR'
 504 !      include 'COMMON.IOUNITS'
 505 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 506 !el      COMMON/SECONDARYS/secstruc
 507       character(len=80) :: line
 508       logical :: errflag
 509 !el      external ilen
 510
 511 !el  local variables
 512       integer :: i,ii,lenpre
 513
 514       allocate(secstruc(nres))
 515
 516 !dr      call getenv_loc('SECPREDFIL',secpred)
 517       lenpre=ilen(prefix)
 518       secpred=prefix(:lenpre)//'.spred'
 519
 520 #if defined(WINIFL) || defined(WINPGI)
 521       open(isecpred,file=secpred,status='old',readonly,shared)
 522 #elif (defined CRAY) || (defined AIX)
 523       open(isecpred,file=secpred,status='old',action='read')
 524 #elif (defined G77)
 525       open(isecpred,file=secpred,status='old')
 526 #else
 527       open(isecpred,file=secpred,status='old',action='read')
 528 #endif
 529 ! read secondary structure prediction from JPRED here!
 530 !      read(isecpred,'(A80)',err=100,end=100) line
 531 !      read(line,'(f10.3)',err=110) ftors
 532        read(isecpred,'(f10.3)',err=110) ftors
 533
 534       write (iout,*) 'FTORS factor =',ftors
 535 ! initialize secstruc to any
 536        do i=1,nres
 537         secstruc(i) ='-'
 538       enddo
 539       ndih_constr=0
 540       ndih_nconstr=0
 541
 542       call read_secstr_pred(isecpred,iout,errflag)
 543       if (errflag) then
 544          write(iout,*)'There is a problem with the list of secondary-',&
 545            'structure prediction'
 546          goto 100
 547       endif
 548 ! 8/13/98 Set limits to generating the dihedral angles
 549       do i=1,nres
 550         phibound(1,i)=-pi
 551         phibound(2,i)=pi
 552       enddo
 553
 554       ii=0
 555       do i=1,nres
 556         if ( secstruc(i) .eq. 'H') then
 557 ! Helix restraints for this residue
 558            ii=ii+1
 559            idih_constr(ii)=i
 560            phi0(ii) = 45.0D0*deg2rad
 561            drange(ii)= 5.0D0*deg2rad
 562            phibound(1,i) = phi0(ii)-drange(ii)
 563            phibound(2,i) = phi0(ii)+drange(ii)
 564         else if (secstruc(i) .eq. 'E') then
 565 ! strand restraints for this residue
 566            ii=ii+1
 567            idih_constr(ii)=i
 568            phi0(ii) = 180.0D0*deg2rad
 569            drange(ii)= 5.0D0*deg2rad
 570            phibound(1,i) = phi0(ii)-drange(ii)
 571            phibound(2,i) = phi0(ii)+drange(ii)
 572         else
 573 ! no restraints for this residue
 574            ndih_nconstr=ndih_nconstr+1
 575            idih_nconstr(ndih_nconstr)=i
 576         endif
 577       enddo
 578       ndih_constr=ii
 579 !      deallocate(secstruc)
 580       return
 581 100   continue
 582       write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
 583 !      deallocate(secstruc)
 584       return
 585 110   continue
 586       write(iout,'(A20)')'Error reading FTORS'
 587 !      deallocate(secstruc)
 588       return
 589       end subroutine secstrp2dihc
 590 !-----------------------------------------------------------------------------
 591       subroutine read_secstr_pred(jin,jout,errors)
 592
 593 !      implicit real*8 (a-h,o-z)
 594 !      INCLUDE 'DIMENSIONS'
 595 !      include 'COMMON.IOUNITS'
 596 !      include 'COMMON.CHAIN'
 597 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 598 !el      COMMON/SECONDARYS/secstruc
 599 !el      EXTERNAL ILEN
 600       character(len=80) :: line,line1   !,ucase
 601       logical :: errflag,errors,blankline
 602
 603 !el  local variables
 604       integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
 605             length_of_chain
 606       errors=.false.
 607       read (jin,'(a)') line
 608       write (jout,'(2a)') '> ',line(1:78)
 609       line1=ucase(line)
 610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
 611 ! correspond to the end-groups.  ADD to the secondary structure prediction "-" for the
 612 ! end-groups in the input file "*.spred"
 613
 614       iseq=1
 615       do while (index(line1,'$END').eq.0)
 616 ! Override commented lines.
 617          ipos=1
 618          blankline=.false.
 619          do while (.not.blankline)
 620             line1=' '
 621             call mykey(line,line1,ipos,blankline,errflag)
 622             if (errflag) write (jout,'(2a)') &
 623        'Error when reading sequence in line: ',line
 624             errors=errors .or. errflag
 625             if (.not. blankline .and. .not. errflag) then
 626                ipos1=2
 627                iend=ilen(line1)
 628 !el               if (iseq.le.maxres) then
 629                   if (line1(1:1).eq.'-' ) then
 630                      secstruc(iseq)=line1(1:1)
 631                   else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
 632                             ( ucase(line1(1:1)).eq.'H' ) ) then
 633                      secstruc(iseq)=ucase(line1(1:1))
 634                   else
 635                      errors=.true.
 636                      write (jout,1010) line1(1:1), iseq
 637                      goto 80
 638                   endif
 639 !el               else
 640 !el                  errors=.true.
 641 !el                  write (jout,1000) iseq,maxres
 642 !el                  goto 80
 643 !el               endif
 644                do while (ipos1.le.iend)
 645
 646                   iseq=iseq+1
 647                   il=1
 648                   ipos1=ipos1+1
 649 !el                  if (iseq.le.maxres) then
 650                      if (line1(ipos1-1:ipos1-1).eq.'-' ) then
 651                         secstruc(iseq)=line1(ipos1-1:ipos1-1)
 652                      else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
 653                            (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
 654                         secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
 655                      else
 656                         errors=.true.
 657                         write (jout,1010) line1(ipos1-1:ipos1-1), iseq
 658                         goto 80
 659                      endif
 660 !el                  else
 661 !el                     errors=.true.
 662 !el                     write (jout,1000) iseq,maxres
 663 !el                     goto 80
 664 !el                  endif
 665                enddo
 666                iseq=iseq+1
 667             endif
 668          enddo
 669          read (jin,'(a)') line
 670          write (jout,'(2a)') '> ',line(1:78)
 671          line1=ucase(line)
 672       enddo
 673
 674 !d    write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
 675
 676 !d check whether the found length of the chain is correct.
 677       length_of_chain=iseq-1
 678       if (length_of_chain .ne. nres) then
 679 !        errors=.true.
 680         write (jout,'(a,i4,a,i4,a)') &
 681        'Error: the number of labels specified in $SEC_STRUC_PRED (' &
 682        ,length_of_chain,') does not match with the number of residues (' &
 683        ,nres,').'
 684       endif
 685    80 continue
 686
 687  1000 format('Error - the number of residues (',i4,&
 688        ') has exceeded maximum (',i4,').')
 689  1010 format ('Error - unrecognized secondary structure label',a4,&
 690        ' in position',i4)
 691       return
 692       end subroutine read_secstr_pred
 693 !#endif
 694 !-----------------------------------------------------------------------------
 695 ! parmread.F
 696 !-----------------------------------------------------------------------------
 697       subroutine parmread
 698
 699       use geometry_data
 700       use energy_data
 701       use control_data, only:maxtor,maxterm
 702       use MD_data
 703       use MPI_data
 704 !el      use map_data
 705       use control, only: getenv_loc
 706 !
 707 ! Read the parameters of the probability distributions of the virtual-bond
 708 ! valence angles and the side chains and energy parameters.
 709 !
 710 ! Important! Energy-term weights ARE NOT read here; they are read from the
 711 ! main input file instead, because NO defaults have yet been set for these
 712 ! parameters.
 713 !
 714 !      implicit real*8 (a-h,o-z)
 715 !      include 'DIMENSIONS'
 716 #ifdef MPI
 717       include "mpif.h"
 718       integer :: IERROR
 719 #endif
 720 !      include 'COMMON.IOUNITS'
 721 !      include 'COMMON.CHAIN'
 722 !      include 'COMMON.INTERACT'
 723 !      include 'COMMON.GEO'
 724 !      include 'COMMON.LOCAL'
 725 !      include 'COMMON.TORSION'
 726 !      include 'COMMON.SCCOR'
 727 !      include 'COMMON.SCROT'
 728 !      include 'COMMON.FFIELD'
 729 !      include 'COMMON.NAMES'
 730 !      include 'COMMON.SBRIDGE'
 731 !      include 'COMMON.MD'
 732 !      include 'COMMON.SETUP'
 733       character(len=1) :: t1,t2,t3
 734       character(len=1) :: onelett(4) = (/"G","A","P","D"/)
 735       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
 736       logical :: lprint,LaTeX
 737       real(kind=8),dimension(3,3,maxlob) :: blower      !(3,3,maxlob)
 738       real(kind=8),dimension(13) :: bN
 739       character(len=3) :: lancuch       !,ucase
 740 !el  local variables
 741       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
 742       integer :: maxinter,junk,kk,ii
 743       real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
 744                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
 745                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
 746                 res1
 747       integer :: ichir1,ichir2
 748 !      real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
 749 !el      allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
 750 !el      allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
 751 !
 752 ! For printing parameters after they are read set the following in the UNRES
 753 ! C-shell script:
 754 !
 755 ! setenv PRINT_PARM YES
 756 !
 757 ! To print parameters in LaTeX format rather than as ASCII tables:
 758 !
 759 ! setenv LATEX YES
 760 !
 761       call getenv_loc("PRINT_PARM",lancuch)
 762       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
 763       call getenv_loc("LATEX",lancuch)
 764       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
 765 !
 766       dwa16=2.0d0**(1.0d0/6.0d0)
 767       itypro=20
 768 ! Assign virtual-bond length
 769       vbl=3.8D0
 770       vblinv=1.0D0/vbl
 771       vblinv2=vblinv*vblinv
 772 !
 773 ! Read the virtual-bond parameters, masses, and moments of inertia
 774 ! and Stokes' radii of the peptide group and side chains
 775 !
 776       allocate(dsc(ntyp1)) !(ntyp1)
 777       allocate(dsc_inv(ntyp1)) !(ntyp1)
 778       allocate(nbondterm(ntyp)) !(ntyp)
 779       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 780       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 781       allocate(msc(ntyp+1)) !(ntyp+1)
 782       allocate(isc(ntyp+1)) !(ntyp+1)
 783       allocate(restok(ntyp+1)) !(ntyp+1)
 784       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 785
 786       dsc(:)=0.0d0
 787       dsc_inv(:)=0.0d0
 788
 789 #ifdef CRYST_BOND
 790       read (ibond,*) vbldp0,akp,mp,ip,pstok
 791       do i=1,ntyp
 792         nbondterm(i)=1
 793         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
 794         dsc(i) = vbldsc0(1,i)
 795         if (i.eq.10) then
 796           dsc_inv(i)=0.0D0
 797         else
 798           dsc_inv(i)=1.0D0/dsc(i)
 799         endif
 800       enddo
 801 #else
 802       read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
 803       do i=1,ntyp
 804         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
 805          j=1,nbondterm(i)),msc(i),isc(i),restok(i)
 806         dsc(i) = vbldsc0(1,i)
 807         if (i.eq.10) then
 808           dsc_inv(i)=0.0D0
 809         else
 810           dsc_inv(i)=1.0D0/dsc(i)
 811         endif
 812       enddo
 813 #endif
 814       if (lprint) then
 815         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
 816         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
 817          'inertia','Pstok'
 818         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
 819         do i=1,ntyp
 820           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
 821             vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
 822           do j=2,nbondterm(i)
 823             write (iout,'(13x,3f10.5)') &
 824               vbldsc0(j,i),aksc(j,i),abond0(j,i)
 825           enddo
 826         enddo
 827       endif
 828 !----------------------------------------------------
 829       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
 830       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))      !(-ntyp:ntyp)
 831       allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
 832       allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
 833       allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
 834       allocate(gthet(3,-ntyp:ntyp))     !(3,-ntyp:ntyp)
 835
 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(-maxthetyp1:maxthetyp1,&
 963         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 964 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 965       allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,&
 966         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 967 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 968 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 969       allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 970         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 971       allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 972         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 973       allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 974         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 975       allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 976         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 977 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
 978 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 979       allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
 980         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 981       allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
 982         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,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) (bN(ii),ii=1,13)
1773         if (lprint) then
1774           write (iout,*) 'Type',i
1775           write (iout,'(a,i2,a,f10.5)') ('bN(',ii,')=',bN(ii),ii=1,13)
1776         endif
1777         B1(1,i)  = bN(3)
1778         B1(2,i)  = bN(5)
1779         B1(1,-i) = bN(3)
1780         B1(2,-i) = -bN(5)
1781 !        b1(1,i)=0.0d0
1782 !        b1(2,i)=0.0d0
1783         B1tilde(1,i) = bN(3)
1784         B1tilde(2,i) =-bN(5)
1785         B1tilde(1,-i) =-bN(3)
1786         B1tilde(2,-i) =bN(5)
1787 !        b1tilde(1,i)=0.0d0
1788 !        b1tilde(2,i)=0.0d0
1789         B2(1,i)  = bN(2)
1790         B2(2,i)  = bN(4)
1791         B2(1,-i)  =bN(2)
1792         B2(2,-i)  =-bN(4)
1793
1794 !        b2(1,i)=0.0d0
1795 !        b2(2,i)=0.0d0
1796         CC(1,1,i)= bN(7)
1797         CC(2,2,i)=-bN(7)
1798         CC(2,1,i)= bN(9)
1799         CC(1,2,i)= bN(9)
1800         CC(1,1,-i)= bN(7)
1801         CC(2,2,-i)=-bN(7)
1802         CC(2,1,-i)=-bN(9)
1803         CC(1,2,-i)=-bN(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)=bN(7)
1809         Ctilde(1,2,i)=bN(9)
1810         Ctilde(2,1,i)=-bN(9)
1811         Ctilde(2,2,i)=bN(7)
1812         Ctilde(1,1,-i)=bN(7)
1813         Ctilde(1,2,-i)=-bN(9)
1814         Ctilde(2,1,-i)=bN(9)
1815         Ctilde(2,2,-i)=bN(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)= bN(6)
1822         DD(2,2,i)=-bN(6)
1823         DD(2,1,i)= bN(8)
1824         DD(1,2,i)= bN(8)
1825         DD(1,1,-i)= bN(6)
1826         DD(2,2,-i)=-bN(6)
1827         DD(2,1,-i)=-bN(8)
1828         DD(1,2,-i)=-bN(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)=bN(6)
1834         Dtilde(1,2,i)=bN(8)
1835         Dtilde(2,1,i)=-bN(8)
1836         Dtilde(2,2,i)=bN(6)
1837         Dtilde(1,1,-i)=bN(6)
1838         Dtilde(1,2,-i)=-bN(8)
1839         Dtilde(2,1,-i)=bN(8)
1840         Dtilde(2,2,-i)=bN(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)= bN(10)+bN(11)
1847         EE(2,2,i)=-bN(10)+bN(11)
1848         EE(2,1,i)= bN(12)-bN(13)
1849         EE(1,2,i)= bN(12)+bN(13)
1850         EE(1,1,-i)= bN(10)+bN(11)
1851         EE(2,2,-i)=-bN(10)+bN(11)
1852         EE(2,1,-i)=-bN(12)+bN(13)
1853         EE(1,2,-i)=-bN(12)-bN(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+1
2329           ishift1=ishift1+1
2330           itype(ires_old)=ntyp1
2331           ibeg=2
2332 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2333           if (unres_pdb) then
2334             do j=1,3
2335               dc(j,ires)=sccor(j,iii)
2336             enddo
2337           else
2338             call sccenter(ires,iii,sccor)
2339           endif
2340           iii=0
2341         endif
2342 ! Read free energy
2343         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2344 ! Fish out the ATOM cards.
2345         if (index(card(1:4),'ATOM').gt.0) then
2346           read (card(12:16),*) atom
2347 !          write (iout,*) "! ",atom," !",ires
2348 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2349           read (card(23:26),*) ires
2350           read (card(18:20),'(a3)') res
2351 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2352 !     &      " ires_old",ires_old
2353 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2354 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2355           if (ires-ishift+ishift1.ne.ires_old) then
2356 ! Calculate the CM of the preceding residue.
2357 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2358             if (ibeg.eq.0) then
2359 !              write (iout,*) "Calculating sidechain center iii",iii
2360               if (unres_pdb) then
2361                 do j=1,3
2362                   dc(j,ires+nres)=sccor(j,iii)
2363                 enddo
2364               else
2365                 call sccenter(ires_old,iii,sccor)
2366               endif
2367               iii=0
2368             endif
2369 ! Start new residue.
2370             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2371               ires=ires_old
2372               cycle
2373             else if (ibeg.eq.1) then
2374               write (iout,*) "BEG ires",ires
2375               ishift=ires-1
2376               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2377                 ishift=ishift-1
2378                 itype(1)=ntyp1
2379               endif
2380               ires=ires-ishift+ishift1
2381               ires_old=ires
2382 !              write (iout,*) "ishift",ishift," ires",ires,&
2383 !               " ires_old",ires_old
2384               ibeg=0
2385             else if (ibeg.eq.2) then
2386 ! Start a new chain
2387               ishift=-ires_old+ires-1 !!!!!
2388               ishift1=ishift1-1    !!!!!
2389 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2390               ires=ires-ishift+ishift1
2391               ires_old=ires
2392               ibeg=0
2393             else
2394               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2395               ires=ires-ishift+ishift1
2396               ires_old=ires
2397             endif
2398             if (res.eq.'ACE' .or. res.eq.'NHE') then
2399               itype(ires)=10
2400             else
2401               itype(ires)=rescode(ires,res,0)
2402             endif
2403           else
2404             ires=ires-ishift+ishift1
2405           endif
2406 !          write (iout,*) "ires_old",ires_old," ires",ires
2407           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2408 !            ishift1=ishift1+1
2409           endif
2410 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2411           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2412              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2413             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2414 !            write (iout,*) "backbone ",atom
2415 #ifdef DEBUG
2416             write (iout,'(2i3,2x,a,3f8.3)') &
2417             ires,itype(ires),res,(c(j,ires),j=1,3)
2418 #endif
2419             iii=iii+1
2420             do j=1,3
2421               sccor(j,iii)=c(j,ires)
2422             enddo
2423 !            write (*,*) card(23:27),ires,itype(ires)
2424           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2425                    atom.ne.'N' .and. atom.ne.'C' .and. &
2426                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2427                    atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2428 !            write (iout,*) "sidechain ",atom
2429             iii=iii+1
2430             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2431           endif
2432         endif
2433       enddo
2434    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2435       if (ires.eq.0) return
2436 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2437 ! system
2438       nres=ires
2439       do i=2,nres-1
2440 !        write (iout,*) i,itype(i)
2441         if (itype(i).eq.ntyp1) then
2442 !          write (iout,*) "dummy",i,itype(i)
2443           do j=1,3
2444             c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2445 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
2446             dc(j,i)=c(j,i)
2447           enddo
2448         endif
2449       enddo
2450 ! Calculate the CM of the last side chain.
2451       if (iii.gt.0)  then
2452       if (unres_pdb) then
2453         do j=1,3
2454           dc(j,ires)=sccor(j,iii)
2455         enddo
2456       else
2457         call sccenter(ires,iii,sccor)
2458       endif
2459       endif
2460 !      nres=ires
2461       nsup=nres
2462       nstart_sup=1
2463       if (itype(nres).ne.10) then
2464         nres=nres+1
2465         itype(nres)=ntyp1
2466         if (unres_pdb) then
2467 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2468           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2469           if (fail) then
2470             e2(1)=0.0d0
2471             e2(2)=1.0d0
2472             e2(3)=0.0d0
2473           endif
2474           do j=1,3
2475             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
2476           enddo
2477         else
2478         do j=1,3
2479           dcj=c(j,nres-2)-c(j,nres-3)
2480           c(j,nres)=c(j,nres-1)+dcj
2481           c(j,2*nres)=c(j,nres)
2482         enddo
2483         endif
2484       endif
2485 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2486 #ifdef WHAM_RUN
2487       if (nres.ne.nres0) then
2488         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2489                        " NRES0=",nres0
2490         stop "Error nres value in WHAM input"
2491       endif
2492 #endif
2493 !---------------------------------
2494 !el reallocate tables
2495 !      do i=1,maxres/3
2496 !       do j=1,2
2497 !         hfrag_alloc(j,i)=hfrag(j,i)
2498 !        enddo
2499 !       do j=1,4
2500 !         bfrag_alloc(j,i)=bfrag(j,i)
2501 !        enddo
2502 !      enddo
2503
2504 !      deallocate(hfrag)
2505 !      deallocate(bfrag)
2506 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
2507 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2508 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2509 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
2510
2511 !      do i=1,nhfrag
2512 !       do j=1,2
2513 !         hfrag(j,i)=hfrag_alloc(j,i)
2514 !        enddo
2515 !      enddo
2516 !      do i=1,nbfrag
2517 !       do j=1,4
2518 !         bfrag(j,i)=bfrag_alloc(j,i)
2519 !        enddo
2520 !      enddo
2521 !el end reallocate tables
2522 !---------------------------------
2523       do i=2,nres-1
2524         do j=1,3
2525           c(j,i+nres)=dc(j,i)
2526         enddo
2527       enddo
2528       do j=1,3
2529         c(j,nres+1)=c(j,1)
2530         c(j,2*nres)=c(j,nres)
2531       enddo
2532       if (itype(1).eq.ntyp1) then
2533         nsup=nsup-1
2534         nstart_sup=2
2535         if (unres_pdb) then
2536 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2537           call refsys(2,3,4,e1,e2,e3,fail)
2538           if (fail) then
2539             e2(1)=0.0d0
2540             e2(2)=1.0d0
2541             e2(3)=0.0d0
2542           endif
2543           do j=1,3
2544             c(j,1)=c(j,2)-3.8d0*e2(j)
2545           enddo
2546         else
2547         do j=1,3
2548           dcj=c(j,4)-c(j,3)
2549           c(j,1)=c(j,2)-dcj
2550           c(j,nres+1)=c(j,1)
2551         enddo
2552         endif
2553       endif
2554 ! Copy the coordinates to reference coordinates
2555 !      do i=1,2*nres
2556 !        do j=1,3
2557 !          cref(j,i)=c(j,i)
2558 !        enddo
2559 !      enddo
2560 ! Calculate internal coordinates.
2561       if (lprn) then
2562       write (iout,'(/a)') &
2563         "Cartesian coordinates of the reference structure"
2564       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2565        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2566       do ires=1,nres
2567         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2568           restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2569           (c(j,ires+nres),j=1,3)
2570       enddo
2571       endif
2572 ! znamy już nres wiec mozna alokowac tablice
2573 ! Calculate internal coordinates.
2574       if(me.eq.king.or..not.out1file)then
2575        write (iout,'(a)') &
2576          "Backbone and SC coordinates as read from the PDB"
2577        do ires=1,nres
2578         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2579           ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2580           (c(j,nres+ires),j=1,3)
2581        enddo
2582       endif
2583
2584       if(.not.allocated(vbld)) then
2585        allocate(vbld(2*nres))
2586        do i=1,2*nres
2587          vbld(i)=0.d0
2588        enddo
2589       endif
2590       if(.not.allocated(vbld_inv)) then
2591        allocate(vbld_inv(2*nres))
2592        do i=1,2*nres
2593          vbld_inv(i)=0.d0
2594        enddo
2595       endif
2596 !!!el
2597       if(.not.allocated(theta)) then
2598         allocate(theta(nres+2))
2599         theta(:)=0.0d0
2600       endif
2601
2602       if(.not.allocated(phi)) allocate(phi(nres+2))
2603       if(.not.allocated(alph)) allocate(alph(nres+2))
2604       if(.not.allocated(omeg)) allocate(omeg(nres+2))
2605       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2606       if(.not.allocated(phiref)) allocate(phiref(nres+2))
2607       if(.not.allocated(costtab)) allocate(costtab(nres))
2608       if(.not.allocated(sinttab)) allocate(sinttab(nres))
2609       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2610       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2611       if(.not.allocated(xxref)) allocate(xxref(nres))
2612       if(.not.allocated(yyref)) allocate(yyref(nres))
2613       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2614       if(.not.allocated(dc_norm)) then
2615 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2616         allocate(dc_norm(3,0:2*nres+2))
2617         dc_norm(:,:)=0.d0
2618       endif
2619
2620       call int_from_cart(.true.,.false.)
2621       call sc_loc_geom(.false.)
2622       do i=1,nres
2623         thetaref(i)=theta(i)
2624         phiref(i)=phi(i)
2625       enddo
2626 !      do i=1,2*nres
2627 !        vbld_inv(i)=0.d0
2628 !        vbld(i)=0.d0
2629 !      enddo
2630
2631       do i=1,nres-1
2632         do j=1,3
2633           dc(j,i)=c(j,i+1)-c(j,i)
2634           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2635         enddo
2636       enddo
2637       do i=2,nres-1
2638         do j=1,3
2639           dc(j,i+nres)=c(j,i+nres)-c(j,i)
2640           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2641         enddo
2642 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2643 !        vbld_inv(i+nres)
2644       enddo
2645 !      call chainbuild
2646 ! Copy the coordinates to reference coordinates
2647 ! Splits to single chain if occurs
2648
2649 !      do i=1,2*nres
2650 !        do j=1,3
2651 !          cref(j,i,cou)=c(j,i)
2652 !        enddo
2653 !      enddo
2654 !
2655       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2656       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2657       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2658 !-----------------------------
2659       kkk=1
2660       lll=0
2661       cou=1
2662       do i=1,nres
2663       lll=lll+1
2664 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2665       if (i.gt.1) then
2666       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2667       chain_length=lll-1
2668       kkk=kkk+1
2669 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2670       lll=1
2671       endif
2672       endif
2673         do j=1,3
2674           cref(j,i,cou)=c(j,i)
2675           cref(j,i+nres,cou)=c(j,i+nres)
2676           if (i.le.nres) then
2677           chain_rep(j,lll,kkk)=c(j,i)
2678           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2679           endif
2680          enddo
2681       enddo
2682       write (iout,*) chain_length
2683       if (chain_length.eq.0) chain_length=nres
2684       do j=1,3
2685       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2686       chain_rep(j,chain_length+nres,symetr) &
2687       =chain_rep(j,chain_length+nres,1)
2688       enddo
2689 ! diagnostic
2690 !       write (iout,*) "spraw lancuchy",chain_length,symetr
2691 !       do i=1,4
2692 !         do kkk=1,chain_length
2693 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2694 !         enddo
2695 !        enddo
2696 ! enddiagnostic
2697 ! makes copy of chains
2698         write (iout,*) "symetr", symetr
2699
2700       if (symetr.gt.1) then
2701        call permut(symetr)
2702        nperm=1
2703        do i=1,symetr
2704        nperm=nperm*i
2705        enddo
2706        do i=1,nperm
2707        write(iout,*) (tabperm(i,kkk),kkk=1,4)
2708        enddo
2709        do i=1,nperm
2710         cou=0
2711         do kkk=1,symetr
2712          icha=tabperm(i,kkk)
2713 !         write (iout,*) i,icha
2714          do lll=1,chain_length
2715           cou=cou+1
2716            if (cou.le.nres) then
2717            do j=1,3
2718             kupa=mod(lll,chain_length)
2719             iprzes=(kkk-1)*chain_length+lll
2720             if (kupa.eq.0) kupa=chain_length
2721 !            write (iout,*) "kupa", kupa
2722             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2723             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2724           enddo
2725           endif
2726          enddo
2727         enddo
2728        enddo
2729        endif
2730 !-koniec robienia kopii
2731 ! diag
2732       do kkk=1,nperm
2733       write (iout,*) "nowa struktura", nperm
2734       do i=1,nres
2735         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2736       cref(2,i,kkk),&
2737       cref(3,i,kkk),cref(1,nres+i,kkk),&
2738       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2739       enddo
2740   100 format (//'              alpha-carbon coordinates       ',&
2741                 '     centroid coordinates'/ &
2742                 '       ', 6X,'X',11X,'Y',11X,'Z', &
2743                                 10X,'X',11X,'Y',11X,'Z')
2744   110 format (a,'(',i3,')',6f12.5)
2745
2746       enddo
2747 !c enddiag
2748       do j=1,nbfrag
2749         do i=1,4
2750          bfrag(i,j)=bfrag(i,j)-ishift
2751         enddo
2752       enddo
2753
2754       do j=1,nhfrag
2755         do i=1,2
2756          hfrag(i,j)=hfrag(i,j)-ishift
2757         enddo
2758       enddo
2759       ishift_pdb=ishift
2760
2761       return
2762       end subroutine readpdb
2763 #if !defined(WHAM_RUN) && !defined(CLUSTER)
2764 !-----------------------------------------------------------------------------
2765 ! readrtns_CSA.F
2766 !-----------------------------------------------------------------------------
2767       subroutine read_control
2768 !
2769 ! Read contorl data
2770 !
2771 !      use geometry_data
2772       use comm_machsw
2773       use energy_data
2774       use control_data
2775       use compare_data
2776       use MCM_data
2777       use map_data
2778       use csa_data
2779       use MD_data
2780       use MPI_data
2781       use random, only: random_init
2782 !      implicit real*8 (a-h,o-z)
2783 !      include 'DIMENSIONS'
2784 #ifdef MP
2785       use prng, only:prng_restart
2786       include 'mpif.h'
2787       logical :: OKRandom!, prng_restart
2788       real(kind=8) :: r1
2789 #endif
2790 !      include 'COMMON.IOUNITS'
2791 !      include 'COMMON.TIME1'
2792 !      include 'COMMON.THREAD'
2793 !      include 'COMMON.SBRIDGE'
2794 !      include 'COMMON.CONTROL'
2795 !      include 'COMMON.MCM'
2796 !      include 'COMMON.MAP'
2797 !      include 'COMMON.HEADER'
2798 !      include 'COMMON.CSA'
2799 !      include 'COMMON.CHAIN'
2800 !      include 'COMMON.MUCA'
2801 !      include 'COMMON.MD'
2802 !      include 'COMMON.FFIELD'
2803 !      include 'COMMON.INTERACT'
2804 !      include 'COMMON.SETUP'
2805 !el      integer :: KDIAG,ICORFL,IXDR
2806 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2807       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2808         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
2809 !      character(len=80) :: ucase
2810       character(len=640) :: controlcard
2811
2812       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2813
2814
2815       nglob_csa=0
2816       eglob_csa=1d99
2817       nmin_csa=0
2818       read (INP,'(a)') titel
2819       call card_concat(controlcard,.true.)
2820 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2821 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2822       call reada(controlcard,'SEED',seed,0.0D0)
2823       call random_init(seed)
2824 ! Set up the time limit (caution! The time must be input in minutes!)
2825       read_cart=index(controlcard,'READ_CART').gt.0
2826       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2827       call readi(controlcard,'SYM',symetr,1)
2828       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2829       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2830       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2831       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2832       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2833       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2834       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2835       call reada(controlcard,'DRMS',drms,0.1D0)
2836       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2837        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2838        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2839        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2840        write (iout,'(a,f10.1)')'DRMS    = ',drms
2841        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2842        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2843       endif
2844       call readi(controlcard,'NZ_START',nz_start,0)
2845       call readi(controlcard,'NZ_END',nz_end,0)
2846       call readi(controlcard,'IZ_SC',iz_sc,0)
2847       timlim=60.0D0*timlim
2848       safety = 60.0d0*safety
2849       timem=timlim
2850       modecalc=0
2851       call reada(controlcard,"T_BATH",t_bath,300.0d0)
2852       minim=(index(controlcard,'MINIMIZE').gt.0)
2853       dccart=(index(controlcard,'CART').gt.0)
2854       overlapsc=(index(controlcard,'OVERLAP').gt.0)
2855       overlapsc=.not.overlapsc
2856       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
2857       searchsc=.not.searchsc
2858       sideadd=(index(controlcard,'SIDEADD').gt.0)
2859       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
2860       outpdb=(index(controlcard,'PDBOUT').gt.0)
2861       outmol2=(index(controlcard,'MOL2OUT').gt.0)
2862       pdbref=(index(controlcard,'PDBREF').gt.0)
2863       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
2864       indpdb=index(controlcard,'PDBSTART')
2865       extconf=(index(controlcard,'EXTCONF').gt.0)
2866       call readi(controlcard,'IPRINT',iprint,0)
2867       call readi(controlcard,'MAXGEN',maxgen,10000)
2868       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
2869       call readi(controlcard,"KDIAG",kdiag,0)
2870       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
2871       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
2872        write (iout,*) "RESCALE_MODE",rescale_mode
2873       split_ene=index(controlcard,'SPLIT_ENE').gt.0
2874       if (index(controlcard,'REGULAR').gt.0.0D0) then
2875         call reada(controlcard,'WEIDIS',weidis,0.1D0)
2876         modecalc=1
2877         refstr=.true.
2878       endif
2879       if (index(controlcard,'CHECKGRAD').gt.0) then
2880         modecalc=5
2881         if (index(controlcard,'CART').gt.0) then
2882           icheckgrad=1
2883         elseif (index(controlcard,'CARINT').gt.0) then
2884           icheckgrad=2
2885         else
2886           icheckgrad=3
2887         endif
2888       elseif (index(controlcard,'THREAD').gt.0) then
2889         modecalc=2
2890         call readi(controlcard,'THREAD',nthread,0)
2891         if (nthread.gt.0) then
2892           call reada(controlcard,'WEIDIS',weidis,0.1D0)
2893         else
2894           if (fg_rank.eq.0) &
2895           write (iout,'(a)')'A number has to follow the THREAD keyword.'
2896           stop 'Error termination in Read_Control.'
2897         endif
2898       else if (index(controlcard,'MCMA').gt.0) then
2899         modecalc=3
2900       else if (index(controlcard,'MCEE').gt.0) then
2901         modecalc=6
2902       else if (index(controlcard,'MULTCONF').gt.0) then
2903         modecalc=4
2904       else if (index(controlcard,'MAP').gt.0) then
2905         modecalc=7
2906         call readi(controlcard,'MAP',nmap,0)
2907       else if (index(controlcard,'CSA').gt.0) then
2908         modecalc=8
2909 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
2910 !rc
2911 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
2912 !rc
2913 !rc        modecalc=9
2914 !fcm      else if (index(controlcard,'MCMF').gt.0) then
2915 !fmc        modecalc=10
2916       else if (index(controlcard,'SOFTREG').gt.0) then
2917         modecalc=11
2918       else if (index(controlcard,'CHECK_BOND').gt.0) then
2919         modecalc=-1
2920       else if (index(controlcard,'TEST').gt.0) then
2921         modecalc=-2
2922       else if (index(controlcard,'MD').gt.0) then
2923         modecalc=12
2924       else if (index(controlcard,'RE ').gt.0) then
2925         modecalc=14
2926       endif
2927
2928       lmuca=index(controlcard,'MUCA').gt.0
2929       call readi(controlcard,'MUCADYN',mucadyn,0)
2930       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
2931       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
2932        then
2933        write (iout,*) 'MUCADYN=',mucadyn
2934        write (iout,*) 'MUCASMOOTH=',muca_smooth
2935       endif
2936
2937       iscode=index(controlcard,'ONE_LETTER')
2938       indphi=index(controlcard,'PHI')
2939       indback=index(controlcard,'BACK')
2940       iranconf=index(controlcard,'RAND_CONF')
2941       i2ndstr=index(controlcard,'USE_SEC_PRED')
2942       gradout=index(controlcard,'GRADOUT').gt.0
2943       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
2944       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
2945       if (me.eq.king .or. .not.out1file ) &
2946         write (iout,*) "DISTCHAINMAX",distchainmax
2947
2948       if(me.eq.king.or..not.out1file) &
2949        write (iout,'(2a)') diagmeth(kdiag),&
2950         ' routine used to diagonalize matrices.'
2951       return
2952       end subroutine read_control
2953 !-----------------------------------------------------------------------------
2954       subroutine read_REMDpar
2955 !
2956 ! Read REMD settings
2957 !
2958 !       use control
2959 !       use energy
2960 !       use geometry
2961       use REMD_data
2962       use MPI_data
2963       use control_data, only:out1file
2964 !      implicit real*8 (a-h,o-z)
2965 !      include 'DIMENSIONS'
2966 !      include 'COMMON.IOUNITS'
2967 !      include 'COMMON.TIME1'
2968 !      include 'COMMON.MD'
2969       use MD_data
2970 !el #ifndef LANG0
2971 !el      include 'COMMON.LANGEVIN'
2972 !el #else
2973 !el      include 'COMMON.LANGEVIN.lang0'
2974 !el #endif
2975 !      include 'COMMON.INTERACT'
2976 !      include 'COMMON.NAMES'
2977 !      include 'COMMON.GEO'
2978 !      include 'COMMON.REMD'
2979 !      include 'COMMON.CONTROL'
2980 !      include 'COMMON.SETUP'
2981 !      character(len=80) :: ucase
2982       character(len=320) :: controlcard
2983       character(len=3200) :: controlcard1
2984       integer :: iremd_m_total
2985 !el local variables
2986       integer :: i
2987 !     real(kind=8) :: var,ene
2988
2989       if(me.eq.king.or..not.out1file) &
2990        write (iout,*) "REMD setup"
2991
2992       call card_concat(controlcard,.true.)
2993       call readi(controlcard,"NREP",nrep,3)
2994       call readi(controlcard,"NSTEX",nstex,1000)
2995       call reada(controlcard,"RETMIN",retmin,10.0d0)
2996       call reada(controlcard,"RETMAX",retmax,1000.0d0)
2997       mremdsync=(index(controlcard,'SYNC').gt.0)
2998       call readi(controlcard,"NSYN",i_sync_step,100)
2999       restart1file=(index(controlcard,'REST1FILE').gt.0)
3000       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3001       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3002       if(max_cache_traj_use.gt.max_cache_traj) &
3003                 max_cache_traj_use=max_cache_traj
3004       if(me.eq.king.or..not.out1file) then
3005 !d       if (traj1file) then
3006 !rc caching is in testing - NTWX is not ignored
3007 !d        write (iout,*) "NTWX value is ignored"
3008 !d        write (iout,*) "  trajectory is stored to one file by master"
3009 !d        write (iout,*) "  before exchange at NSTEX intervals"
3010 !d       endif
3011        write (iout,*) "NREP= ",nrep
3012        write (iout,*) "NSTEX= ",nstex
3013        write (iout,*) "SYNC= ",mremdsync
3014        write (iout,*) "NSYN= ",i_sync_step
3015        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3016       endif
3017       remd_tlist=.false.
3018       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3019       if (index(controlcard,'TLIST').gt.0) then
3020          remd_tlist=.true.
3021          call card_concat(controlcard1,.true.)
3022          read(controlcard1,*) (remd_t(i),i=1,nrep)
3023          if(me.eq.king.or..not.out1file) &
3024           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3025       endif
3026       remd_mlist=.false.
3027       if (index(controlcard,'MLIST').gt.0) then
3028          remd_mlist=.true.
3029          call card_concat(controlcard1,.true.)
3030          read(controlcard1,*) (remd_m(i),i=1,nrep)
3031          if(me.eq.king.or..not.out1file) then
3032           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3033           iremd_m_total=0
3034           do i=1,nrep
3035            iremd_m_total=iremd_m_total+remd_m(i)
3036           enddo
3037           write (iout,*) 'Total number of replicas ',iremd_m_total
3038          endif
3039       endif
3040       if(me.eq.king.or..not.out1file) &
3041        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3042       return
3043       end subroutine read_REMDpar
3044 !-----------------------------------------------------------------------------
3045       subroutine read_MDpar
3046 !
3047 ! Read MD settings
3048 !
3049       use control_data, only: r_cut,rlamb,out1file
3050       use energy_data
3051       use geometry_data, only: pi
3052       use MPI_data
3053 !      implicit real*8 (a-h,o-z)
3054 !      include 'DIMENSIONS'
3055 !      include 'COMMON.IOUNITS'
3056 !      include 'COMMON.TIME1'
3057 !      include 'COMMON.MD'
3058       use MD_data
3059 !el #ifndef LANG0
3060 !el      include 'COMMON.LANGEVIN'
3061 !el #else
3062 !el      include 'COMMON.LANGEVIN.lang0'
3063 !el #endif
3064 !      include 'COMMON.INTERACT'
3065 !      include 'COMMON.NAMES'
3066 !      include 'COMMON.GEO'
3067 !      include 'COMMON.SETUP'
3068 !      include 'COMMON.CONTROL'
3069 !      include 'COMMON.SPLITELE'
3070 !      character(len=80) :: ucase
3071       character(len=320) :: controlcard
3072 !el local variables
3073       integer :: i
3074       real(kind=8) :: eta
3075
3076       call card_concat(controlcard,.true.)
3077       call readi(controlcard,"NSTEP",n_timestep,1000000)
3078       call readi(controlcard,"NTWE",ntwe,100)
3079       call readi(controlcard,"NTWX",ntwx,1000)
3080       call reada(controlcard,"DT",d_time,1.0d-1)
3081       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3082       call reada(controlcard,"DAMAX",damax,1.0d1)
3083       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3084       call readi(controlcard,"LANG",lang,0)
3085       RESPA = index(controlcard,"RESPA") .gt. 0
3086       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3087       ntime_split0=ntime_split
3088       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3089       ntime_split0=ntime_split
3090       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3091       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3092       rest = index(controlcard,"REST").gt.0
3093       tbf = index(controlcard,"TBF").gt.0
3094       usampl = index(controlcard,"USAMPL").gt.0
3095       mdpdb = index(controlcard,"MDPDB").gt.0
3096       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3097       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3098       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3099       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3100       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3101       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3102       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3103       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3104       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3105       large = index(controlcard,"LARGE").gt.0
3106       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3107       rattle = index(controlcard,"RATTLE").gt.0
3108 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3109       nset=0
3110       if(usampl) then
3111         call read_fragments
3112       endif
3113
3114       if(me.eq.king.or..not.out1file) then
3115        write (iout,*)
3116        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3117        write (iout,*)
3118        write (iout,'(a)') "The units are:"
3119        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3120        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3121         " acceleration: angstrom/(48.9 fs)**2"
3122        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3123        write (iout,*)
3124        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3125        write (iout,'(a60,f10.5,a)') &
3126         "Initial time step of numerical integration:",d_time,&
3127         " natural units"
3128        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3129        if (RESPA) then
3130         write (iout,'(2a,i4,a)') &
3131           "A-MTS algorithm used; initial time step for fast-varying",&
3132           " short-range forces split into",ntime_split," steps."
3133         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3134          r_cut," lambda",rlamb
3135        endif
3136        write (iout,'(2a,f10.5)') &
3137         "Maximum acceleration threshold to reduce the time step",&
3138         "/increase split number:",damax
3139        write (iout,'(2a,f10.5)') &
3140         "Maximum predicted energy drift to reduce the timestep",&
3141         "/increase split number:",edriftmax
3142        write (iout,'(a60,f10.5)') &
3143        "Maximum velocity threshold to reduce velocities:",dvmax
3144        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3145        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3146        if (rattle) write (iout,'(a60)') &
3147         "Rattle algorithm used to constrain the virtual bonds"
3148       endif
3149       reset_fricmat=1000
3150       if (lang.gt.0) then
3151         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3152         call reada(controlcard,"RWAT",rwat,1.4d0)
3153         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3154         surfarea=index(controlcard,"SURFAREA").gt.0
3155         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3156         if(me.eq.king.or..not.out1file)then
3157          write (iout,'(/a,$)') "Langevin dynamics calculation"
3158          if (lang.eq.1) then
3159           write (iout,'(a/)') &
3160             " with direct integration of Langevin equations"
3161          else if (lang.eq.2) then
3162           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3163          else if (lang.eq.3) then
3164           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3165          else if (lang.eq.4) then
3166           write (iout,'(a/)') " in overdamped mode"
3167          else
3168           write (iout,'(//a,i5)') &
3169             "=========== ERROR: Unknown Langevin dynamics mode:",lang
3170           stop
3171          endif
3172          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3173          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3174          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3175          write (iout,'(a60,f10.5)') &
3176          "Scaling factor of the friction forces:",scal_fric
3177          if (surfarea) write (iout,'(2a,i10,a)') &
3178            "Friction coefficients will be scaled by solvent-accessible",&
3179            " surface area every",reset_fricmat," steps."
3180         endif
3181 ! Calculate friction coefficients and bounds of stochastic forces
3182         eta=6*pi*cPoise*etawat
3183         if(me.eq.king.or..not.out1file) &
3184          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3185           eta
3186         gamp=scal_fric*(pstok+rwat)*eta
3187         stdfp=dsqrt(2*Rb*t_bath/d_time)
3188         allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3189         do i=1,ntyp
3190           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3191           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3192         enddo
3193         if(me.eq.king.or..not.out1file)then
3194          write (iout,'(/2a/)') &
3195          "Radii of site types and friction coefficients and std's of",&
3196          " stochastic forces of fully exposed sites"
3197          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3198          do i=1,ntyp
3199           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3200            gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3201          enddo
3202         endif
3203       else if (tbf) then
3204         if(me.eq.king.or..not.out1file)then
3205          write (iout,'(a)') "Berendsen bath calculation"
3206          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3207          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3208          if (reset_moment) &
3209          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3210          count_reset_moment," steps"
3211          if (reset_vel) &
3212           write (iout,'(a,i10,a)') &
3213           "Velocities will be reset at random every",count_reset_vel,&
3214          " steps"
3215         endif
3216       else
3217         if(me.eq.king.or..not.out1file) &
3218          write (iout,'(a31)') "Microcanonical mode calculation"
3219       endif
3220       if(me.eq.king.or..not.out1file)then
3221        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3222        if (usampl) then
3223           write(iout,*) "MD running with constraints."
3224           write(iout,*) "Equilibration time ", eq_time, " mtus."
3225           write(iout,*) "Constraining ", nfrag," fragments."
3226           write(iout,*) "Length of each fragment, weight and q0:"
3227           do iset=1,nset
3228            write (iout,*) "Set of restraints #",iset
3229            do i=1,nfrag
3230               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3231                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3232            enddo
3233            write(iout,*) "constraints between ", npair, "fragments."
3234            write(iout,*) "constraint pairs, weights and q0:"
3235            do i=1,npair
3236             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3237                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3238            enddo
3239            write(iout,*) "angle constraints within ", nfrag_back,&
3240             "backbone fragments."
3241            write(iout,*) "fragment, weights:"
3242            do i=1,nfrag_back
3243             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3244                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3245                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3246            enddo
3247           enddo
3248         iset=mod(kolor,nset)+1
3249        endif
3250       endif
3251       if(me.eq.king.or..not.out1file) &
3252        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3253       return
3254       end subroutine read_MDpar
3255 !-----------------------------------------------------------------------------
3256       subroutine map_read
3257
3258       use map_data
3259 !      implicit real*8 (a-h,o-z)
3260 !      include 'DIMENSIONS'
3261 !      include 'COMMON.MAP'
3262 !      include 'COMMON.IOUNITS'
3263       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3264       character(len=80) :: mapcard      !,ucase
3265 !el local variables
3266       integer :: imap
3267 !     real(kind=8) :: var,ene
3268
3269       do imap=1,nmap
3270         read (inp,'(a)') mapcard
3271         mapcard=ucase(mapcard)
3272         if (index(mapcard,'PHI').gt.0) then
3273           kang(imap)=1
3274         else if (index(mapcard,'THE').gt.0) then
3275           kang(imap)=2
3276         else if (index(mapcard,'ALP').gt.0) then
3277           kang(imap)=3
3278         else if (index(mapcard,'OME').gt.0) then
3279           kang(imap)=4
3280         else
3281           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3282           stop 'Error - illegal variable spec in MAP card.'
3283         endif
3284         call readi (mapcard,'RES1',res1(imap),0)
3285         call readi (mapcard,'RES2',res2(imap),0)
3286         if (res1(imap).eq.0) then
3287           res1(imap)=res2(imap)
3288         else if (res2(imap).eq.0) then
3289           res2(imap)=res1(imap)
3290         endif
3291         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3292           write (iout,'(a)') &
3293           'Error - illegal definition of variable group in MAP.'
3294           stop 'Error - illegal definition of variable group in MAP.'
3295         endif
3296         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3297         call reada(mapcard,'TO',ang_to(imap),0.0D0)
3298         call readi(mapcard,'NSTEP',nstep(imap),0)
3299         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3300           write (iout,'(a)') &
3301            'Illegal boundary and/or step size specification in MAP.'
3302           stop 'Illegal boundary and/or step size specification in MAP.'
3303         endif
3304       enddo ! imap
3305       return
3306       end subroutine map_read
3307 !-----------------------------------------------------------------------------
3308       subroutine csaread
3309
3310       use control_data, only: vdisulf
3311       use csa_data
3312 !      implicit real*8 (a-h,o-z)
3313 !      include 'DIMENSIONS'
3314 !      include 'COMMON.IOUNITS'
3315 !      include 'COMMON.GEO'
3316 !      include 'COMMON.CSA'
3317 !      include 'COMMON.BANK'
3318 !      include 'COMMON.CONTROL'
3319 !      character(len=80) :: ucase
3320       character(len=620) :: mcmcard
3321 !el local variables
3322 !     integer :: ntf,ik,iw_pdb
3323 !     real(kind=8) :: var,ene
3324
3325       call card_concat(mcmcard,.true.)
3326
3327       call readi(mcmcard,'NCONF',nconf,50)
3328       call readi(mcmcard,'NADD',nadd,0)
3329       call readi(mcmcard,'JSTART',jstart,1)
3330       call readi(mcmcard,'JEND',jend,1)
3331       call readi(mcmcard,'NSTMAX',nstmax,500000)
3332       call readi(mcmcard,'N0',n0,1)
3333       call readi(mcmcard,'N1',n1,6)
3334       call readi(mcmcard,'N2',n2,4)
3335       call readi(mcmcard,'N3',n3,0)
3336       call readi(mcmcard,'N4',n4,0)
3337       call readi(mcmcard,'N5',n5,0)
3338       call readi(mcmcard,'N6',n6,10)
3339       call readi(mcmcard,'N7',n7,0)
3340       call readi(mcmcard,'N8',n8,0)
3341       call readi(mcmcard,'N9',n9,0)
3342       call readi(mcmcard,'N14',n14,0)
3343       call readi(mcmcard,'N15',n15,0)
3344       call readi(mcmcard,'N16',n16,0)
3345       call readi(mcmcard,'N17',n17,0)
3346       call readi(mcmcard,'N18',n18,0)
3347
3348       vdisulf=(index(mcmcard,'DYNSS').gt.0)
3349
3350       call readi(mcmcard,'NDIFF',ndiff,2)
3351       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3352       call readi(mcmcard,'IS1',is1,1)
3353       call readi(mcmcard,'IS2',is2,8)
3354       call readi(mcmcard,'NRAN0',nran0,4)
3355       call readi(mcmcard,'NRAN1',nran1,2)
3356       call readi(mcmcard,'IRR',irr,1)
3357       call readi(mcmcard,'NSEED',nseed,20)
3358       call readi(mcmcard,'NTOTAL',ntotal,10000)
3359       call reada(mcmcard,'CUT1',cut1,2.0d0)
3360       call reada(mcmcard,'CUT2',cut2,5.0d0)
3361       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3362       call readi(mcmcard,'ICMAX',icmax,3)
3363       call readi(mcmcard,'IRESTART',irestart,0)
3364 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
3365       ntbankm=0
3366 !!bankt
3367       call reada(mcmcard,'DELE',dele,20.0d0)
3368       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3369       call readi(mcmcard,'IREF',iref,0)
3370       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3371       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3372       call readi(mcmcard,'NCONF_IN',nconf_in,0)
3373       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3374       write (iout,*) "NCONF_IN",nconf_in
3375       return
3376       end subroutine csaread
3377 !-----------------------------------------------------------------------------
3378       subroutine mcmread
3379
3380       use mcm_data
3381       use control_data, only: MaxMoveType
3382       use MD_data
3383       use minim_data
3384 !      implicit real*8 (a-h,o-z)
3385 !      include 'DIMENSIONS'
3386 !      include 'COMMON.MCM'
3387 !      include 'COMMON.MCE'
3388 !      include 'COMMON.IOUNITS'
3389 !      character(len=80) :: ucase
3390       character(len=320) :: mcmcard
3391 !el local variables
3392       integer :: i
3393 !     real(kind=8) :: var,ene
3394
3395       call card_concat(mcmcard,.true.)
3396       call readi(mcmcard,'MAXACC',maxacc,100)
3397       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3398       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3399       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3400       call readi(mcmcard,'MAXREPM',maxrepm,200)
3401       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3402       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3403       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3404       call reada(mcmcard,'E_UP',e_up,5.0D0)
3405       call reada(mcmcard,'DELTE',delte,0.1D0)
3406       call readi(mcmcard,'NSWEEP',nsweep,5)
3407       call readi(mcmcard,'NSTEPH',nsteph,0)
3408       call readi(mcmcard,'NSTEPC',nstepc,0)
3409       call reada(mcmcard,'TMIN',tmin,298.0D0)
3410       call reada(mcmcard,'TMAX',tmax,298.0D0)
3411       call readi(mcmcard,'NWINDOW',nwindow,0)
3412       call readi(mcmcard,'PRINT_MC',print_mc,0)
3413       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3414       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3415       ent_read=(index(mcmcard,'ENT_READ').gt.0)
3416       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3417       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3418       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3419       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3420       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3421       if (nwindow.gt.0) then
3422         allocate(winstart(nwindow))     !!el (maxres)
3423         allocate(winend(nwindow))       !!el
3424         allocate(winlen(nwindow))       !!el
3425         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3426         do i=1,nwindow
3427           winlen(i)=winend(i)-winstart(i)+1
3428         enddo
3429       endif
3430       if (tmax.lt.tmin) tmax=tmin
3431       if (tmax.eq.tmin) then
3432         nstepc=0
3433         nsteph=0
3434       endif
3435       if (nstepc.gt.0 .and. nsteph.gt.0) then
3436         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3437         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3438       endif
3439       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3440 ! Probabilities of different move types
3441       sumpro_type(0)=0.0D0
3442       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3443       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3444       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3445       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
3446       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3447       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3448       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3449       do i=1,MaxMoveType
3450         print *,'i',i,' sumprotype',sumpro_type(i)
3451         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3452         print *,'i',i,' sumprotype',sumpro_type(i)
3453       enddo
3454       return
3455       end subroutine mcmread
3456 !-----------------------------------------------------------------------------
3457       subroutine read_minim
3458
3459       use minim_data
3460 !      implicit real*8 (a-h,o-z)
3461 !      include 'DIMENSIONS'
3462 !      include 'COMMON.MINIM'
3463 !      include 'COMMON.IOUNITS'
3464 !      character(len=80) :: ucase
3465       character(len=320) :: minimcard
3466 !el local variables
3467 !     integer :: ntf,ik,iw_pdb
3468 !     real(kind=8) :: var,ene
3469
3470       call card_concat(minimcard,.true.)
3471       call readi(minimcard,'MAXMIN',maxmin,2000)
3472       call readi(minimcard,'MAXFUN',maxfun,5000)
3473       call readi(minimcard,'MINMIN',minmin,maxmin)
3474       call readi(minimcard,'MINFUN',minfun,maxmin)
3475       call reada(minimcard,'TOLF',tolf,1.0D-2)
3476       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3477       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3478       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3479       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3480       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3481                'Options in energy minimization:'
3482       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3483        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3484        'MinMin:',MinMin,' MinFun:',MinFun,&
3485        ' TolF:',TolF,' RTolF:',RTolF
3486       return
3487       end subroutine read_minim
3488 !-----------------------------------------------------------------------------
3489       subroutine openunits
3490
3491       use energy_data, only: usampl
3492       use csa_data
3493       use MPI_data
3494       use control_data, only:out1file
3495       use control, only: getenv_loc
3496 !      implicit real*8 (a-h,o-z)
3497 !      include 'DIMENSIONS'
3498 #ifdef MPI
3499       include 'mpif.h'
3500       character(len=16) :: form,nodename
3501       integer :: nodelen,ierror,npos
3502 #endif
3503 !      include 'COMMON.SETUP'
3504 !      include 'COMMON.IOUNITS'
3505 !      include 'COMMON.MD'
3506 !      include 'COMMON.CONTROL'
3507       integer :: lenpre,lenpot,lentmp   !,ilen
3508 !el      external ilen
3509       character(len=3) :: out1file_text !,ucase
3510       character(len=3) :: ll
3511 !el      external ucase
3512 !el local variables
3513 !     integer :: ntf,ik,iw_pdb
3514 !     real(kind=8) :: var,ene
3515 !
3516 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3517       call getenv_loc("PREFIX",prefix)
3518       pref_orig = prefix
3519       call getenv_loc("POT",pot)
3520       call getenv_loc("DIRTMP",tmpdir)
3521       call getenv_loc("CURDIR",curdir)
3522       call getenv_loc("OUT1FILE",out1file_text)
3523 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3524       out1file_text=ucase(out1file_text)
3525       if (out1file_text(1:1).eq."Y") then
3526         out1file=.true.
3527       else
3528         out1file=fg_rank.gt.0
3529       endif
3530       lenpre=ilen(prefix)
3531       lenpot=ilen(pot)
3532       lentmp=ilen(tmpdir)
3533       if (lentmp.gt.0) then
3534           write (*,'(80(1h!))')
3535           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
3536           write (*,'(80(1h!))')
3537           write (*,*)"All output files will be on node /tmp directory."
3538 #ifdef MPI
3539         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3540         if (me.eq.king) then
3541           write (*,*) "The master node is ",nodename
3542         else if (fg_rank.eq.0) then
3543           write (*,*) "I am the CG slave node ",nodename
3544         else
3545           write (*,*) "I am the FG slave node ",nodename
3546         endif
3547 #endif
3548         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3549         lenpre = lentmp+lenpre+1
3550       endif
3551       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3552 ! Get the names and open the input files
3553 #if defined(WINIFL) || defined(WINPGI)
3554       open(1,file=pref_orig(:ilen(pref_orig))// &
3555         '.inp',status='old',readonly,shared)
3556        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3557 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3558 ! Get parameter filenames and open the parameter files.
3559       call getenv_loc('BONDPAR',bondname)
3560       open (ibond,file=bondname,status='old',readonly,shared)
3561       call getenv_loc('THETPAR',thetname)
3562       open (ithep,file=thetname,status='old',readonly,shared)
3563       call getenv_loc('ROTPAR',rotname)
3564       open (irotam,file=rotname,status='old',readonly,shared)
3565       call getenv_loc('TORPAR',torname)
3566       open (itorp,file=torname,status='old',readonly,shared)
3567       call getenv_loc('TORDPAR',tordname)
3568       open (itordp,file=tordname,status='old',readonly,shared)
3569       call getenv_loc('FOURIER',fouriername)
3570       open (ifourier,file=fouriername,status='old',readonly,shared)
3571       call getenv_loc('ELEPAR',elename)
3572       open (ielep,file=elename,status='old',readonly,shared)
3573       call getenv_loc('SIDEPAR',sidename)
3574       open (isidep,file=sidename,status='old',readonly,shared)
3575 #elif (defined CRAY) || (defined AIX)
3576       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3577         action='read')
3578 !      print *,"Processor",myrank," opened file 1"
3579       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3580 !      print *,"Processor",myrank," opened file 9"
3581 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3582 ! Get parameter filenames and open the parameter files.
3583       call getenv_loc('BONDPAR',bondname)
3584       open (ibond,file=bondname,status='old',action='read')
3585 !      print *,"Processor",myrank," opened file IBOND"
3586       call getenv_loc('THETPAR',thetname)
3587       open (ithep,file=thetname,status='old',action='read')
3588 !      print *,"Processor",myrank," opened file ITHEP"
3589       call getenv_loc('ROTPAR',rotname)
3590       open (irotam,file=rotname,status='old',action='read')
3591 !      print *,"Processor",myrank," opened file IROTAM"
3592       call getenv_loc('TORPAR',torname)
3593       open (itorp,file=torname,status='old',action='read')
3594 !      print *,"Processor",myrank," opened file ITORP"
3595       call getenv_loc('TORDPAR',tordname)
3596       open (itordp,file=tordname,status='old',action='read')
3597 !      print *,"Processor",myrank," opened file ITORDP"
3598       call getenv_loc('SCCORPAR',sccorname)
3599       open (isccor,file=sccorname,status='old',action='read')
3600 !      print *,"Processor",myrank," opened file ISCCOR"
3601       call getenv_loc('FOURIER',fouriername)
3602       open (ifourier,file=fouriername,status='old',action='read')
3603 !      print *,"Processor",myrank," opened file IFOURIER"
3604       call getenv_loc('ELEPAR',elename)
3605       open (ielep,file=elename,status='old',action='read')
3606 !      print *,"Processor",myrank," opened file IELEP"
3607       call getenv_loc('SIDEPAR',sidename)
3608       open (isidep,file=sidename,status='old',action='read')
3609 !      print *,"Processor",myrank," opened file ISIDEP"
3610 !      print *,"Processor",myrank," opened parameter files"
3611 #elif (defined G77)
3612       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3613       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3614 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3615 ! Get parameter filenames and open the parameter files.
3616       call getenv_loc('BONDPAR',bondname)
3617       open (ibond,file=bondname,status='old')
3618       call getenv_loc('THETPAR',thetname)
3619       open (ithep,file=thetname,status='old')
3620       call getenv_loc('ROTPAR',rotname)
3621       open (irotam,file=rotname,status='old')
3622       call getenv_loc('TORPAR',torname)
3623       open (itorp,file=torname,status='old')
3624       call getenv_loc('TORDPAR',tordname)
3625       open (itordp,file=tordname,status='old')
3626       call getenv_loc('SCCORPAR',sccorname)
3627       open (isccor,file=sccorname,status='old')
3628       call getenv_loc('FOURIER',fouriername)
3629       open (ifourier,file=fouriername,status='old')
3630       call getenv_loc('ELEPAR',elename)
3631       open (ielep,file=elename,status='old')
3632       call getenv_loc('SIDEPAR',sidename)
3633       open (isidep,file=sidename,status='old')
3634 #else
3635       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3636         readonly)
3637        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3638 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3639 ! Get parameter filenames and open the parameter files.
3640       call getenv_loc('BONDPAR',bondname)
3641       open (ibond,file=bondname,status='old',action='read')
3642       call getenv_loc('THETPAR',thetname)
3643       open (ithep,file=thetname,status='old',action='read')
3644       call getenv_loc('ROTPAR',rotname)
3645       open (irotam,file=rotname,status='old',action='read')
3646       call getenv_loc('TORPAR',torname)
3647       open (itorp,file=torname,status='old',action='read')
3648       call getenv_loc('TORDPAR',tordname)
3649       open (itordp,file=tordname,status='old',action='read')
3650       call getenv_loc('SCCORPAR',sccorname)
3651       open (isccor,file=sccorname,status='old',action='read')
3652 #ifndef CRYST_THETA
3653       call getenv_loc('THETPARPDB',thetname_pdb)
3654       print *,"thetname_pdb ",thetname_pdb
3655       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3656       print *,ithep_pdb," opened"
3657 #endif
3658       call getenv_loc('FOURIER',fouriername)
3659       open (ifourier,file=fouriername,status='old',readonly)
3660       call getenv_loc('ELEPAR',elename)
3661       open (ielep,file=elename,status='old',readonly)
3662       call getenv_loc('SIDEPAR',sidename)
3663       open (isidep,file=sidename,status='old',readonly)
3664 #ifndef CRYST_SC
3665       call getenv_loc('ROTPARPDB',rotname_pdb)
3666       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3667 #endif
3668 #endif
3669 #ifndef OLDSCP
3670 !
3671 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3672 ! Use -DOLDSCP to use hard-coded constants instead.
3673 !
3674       call getenv_loc('SCPPAR',scpname)
3675 #if defined(WINIFL) || defined(WINPGI)
3676       open (iscpp,file=scpname,status='old',readonly,shared)
3677 #elif (defined CRAY)  || (defined AIX)
3678       open (iscpp,file=scpname,status='old',action='read')
3679 #elif (defined G77)
3680       open (iscpp,file=scpname,status='old')
3681 #else
3682       open (iscpp,file=scpname,status='old',action='read')
3683 #endif
3684 #endif
3685       call getenv_loc('PATTERN',patname)
3686 #if defined(WINIFL) || defined(WINPGI)
3687       open (icbase,file=patname,status='old',readonly,shared)
3688 #elif (defined CRAY)  || (defined AIX)
3689       open (icbase,file=patname,status='old',action='read')
3690 #elif (defined G77)
3691       open (icbase,file=patname,status='old')
3692 #else
3693       open (icbase,file=patname,status='old',action='read')
3694 #endif
3695 #ifdef MPI
3696 ! Open output file only for CG processes
3697 !      print *,"Processor",myrank," fg_rank",fg_rank
3698       if (fg_rank.eq.0) then
3699
3700       if (nodes.eq.1) then
3701         npos=3
3702       else
3703         npos = dlog10(dfloat(nodes-1))+1
3704       endif
3705       if (npos.lt.3) npos=3
3706       write (liczba,'(i1)') npos
3707       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3708         //')'
3709       write (liczba,form) me
3710       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3711         liczba(:ilen(liczba))
3712       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3713         //'.int'
3714       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3715         //'.pdb'
3716       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3717         liczba(:ilen(liczba))//'.mol2'
3718       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3719         liczba(:ilen(liczba))//'.stat'
3720       if (lentmp.gt.0) &
3721         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3722             //liczba(:ilen(liczba))//'.stat')
3723       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3724         //'.rst'
3725       if(usampl) then
3726           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3727        liczba(:ilen(liczba))//'.const'
3728       endif
3729
3730       endif
3731 #else
3732       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3733       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3734       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3735       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3736       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3737       if (lentmp.gt.0) &
3738         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3739           '.stat')
3740       rest2name=prefix(:ilen(prefix))//'.rst'
3741       if(usampl) then
3742          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3743       endif
3744 #endif
3745 #if defined(AIX) || defined(PGI)
3746       if (me.eq.king .or. .not. out1file) &
3747          open(iout,file=outname,status='unknown')
3748 #ifdef DEBUG
3749       if (fg_rank.gt.0) then
3750         write (liczba,'(i3.3)') myrank/nfgtasks
3751         write (ll,'(bz,i3.3)') fg_rank
3752         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3753          status='unknown')
3754       endif
3755 #endif
3756       if(me.eq.king) then
3757        open(igeom,file=intname,status='unknown',position='append')
3758        open(ipdb,file=pdbname,status='unknown')
3759        open(imol2,file=mol2name,status='unknown')
3760        open(istat,file=statname,status='unknown',position='append')
3761       else
3762 !1out       open(iout,file=outname,status='unknown')
3763       endif
3764 #else
3765       if (me.eq.king .or. .not.out1file) &
3766           open(iout,file=outname,status='unknown')
3767 #ifdef DEBUG
3768       if (fg_rank.gt.0) then
3769         write (liczba,'(i3.3)') myrank/nfgtasks
3770         write (ll,'(bz,i3.3)') fg_rank
3771         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3772          status='unknown')
3773       endif
3774 #endif
3775       if(me.eq.king) then
3776        open(igeom,file=intname,status='unknown',access='append')
3777        open(ipdb,file=pdbname,status='unknown')
3778        open(imol2,file=mol2name,status='unknown')
3779        open(istat,file=statname,status='unknown',access='append')
3780       else
3781 !1out       open(iout,file=outname,status='unknown')
3782       endif
3783 #endif
3784       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3785       csa_seed=prefix(:lenpre)//'.CSA.seed'
3786       csa_history=prefix(:lenpre)//'.CSA.history'
3787       csa_bank=prefix(:lenpre)//'.CSA.bank'
3788       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3789       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3790       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3791 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3792       csa_int=prefix(:lenpre)//'.int'
3793       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3794       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3795       csa_in=prefix(:lenpre)//'.CSA.in'
3796 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3797 ! Write file names
3798       if (me.eq.king)then
3799       write (iout,'(80(1h-))')
3800       write (iout,'(30x,a)') "FILE ASSIGNMENT"
3801       write (iout,'(80(1h-))')
3802       write (iout,*) "Input file                      : ",&
3803         pref_orig(:ilen(pref_orig))//'.inp'
3804       write (iout,*) "Output file                     : ",&
3805         outname(:ilen(outname))
3806       write (iout,*)
3807       write (iout,*) "Sidechain potential file        : ",&
3808         sidename(:ilen(sidename))
3809 #ifndef OLDSCP
3810       write (iout,*) "SCp potential file              : ",&
3811         scpname(:ilen(scpname))
3812 #endif
3813       write (iout,*) "Electrostatic potential file    : ",&
3814         elename(:ilen(elename))
3815       write (iout,*) "Cumulant coefficient file       : ",&
3816         fouriername(:ilen(fouriername))
3817       write (iout,*) "Torsional parameter file        : ",&
3818         torname(:ilen(torname))
3819       write (iout,*) "Double torsional parameter file : ",&
3820         tordname(:ilen(tordname))
3821       write (iout,*) "SCCOR parameter file : ",&
3822         sccorname(:ilen(sccorname))
3823       write (iout,*) "Bond & inertia constant file    : ",&
3824         bondname(:ilen(bondname))
3825       write (iout,*) "Bending parameter file          : ",&
3826         thetname(:ilen(thetname))
3827       write (iout,*) "Rotamer parameter file          : ",&
3828         rotname(:ilen(rotname))
3829 !el----
3830 #ifndef CRYST_THETA
3831       write (iout,*) "Thetpdb parameter file          : ",&
3832         thetname_pdb(:ilen(thetname_pdb))
3833 #endif
3834 !el
3835       write (iout,*) "Threading database              : ",&
3836         patname(:ilen(patname))
3837       if (lentmp.ne.0) &
3838       write (iout,*)" DIRTMP                          : ",&
3839         tmpdir(:lentmp)
3840       write (iout,'(80(1h-))')
3841       endif
3842       return
3843       end subroutine openunits
3844 !-----------------------------------------------------------------------------
3845       subroutine readrst
3846
3847       use geometry_data, only: nres,dc
3848       use energy_data, only: usampl,iset
3849       use MD_data
3850 !      implicit real*8 (a-h,o-z)
3851 !      include 'DIMENSIONS'
3852 !      include 'COMMON.CHAIN'
3853 !      include 'COMMON.IOUNITS'
3854 !      include 'COMMON.MD'
3855 !el local variables
3856       integer ::i,j
3857 !     real(kind=8) :: var,ene
3858
3859       open(irest2,file=rest2name,status='unknown')
3860       read(irest2,*) totT,EK,potE,totE,t_bath
3861       do i=1,2*nres
3862          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
3863       enddo
3864       do i=1,2*nres
3865          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
3866       enddo
3867       if(usampl) then
3868              read (irest2,*) iset
3869       endif
3870       close(irest2)
3871       return
3872       end subroutine readrst
3873 !-----------------------------------------------------------------------------
3874       subroutine read_fragments
3875
3876       use energy_data
3877 !      use geometry
3878       use control_data, only:out1file
3879       use MD_data
3880       use MPI_data
3881 !      implicit real*8 (a-h,o-z)
3882 !      include 'DIMENSIONS'
3883 #ifdef MPI
3884       include 'mpif.h'
3885 #endif
3886 !      include 'COMMON.SETUP'
3887 !      include 'COMMON.CHAIN'
3888 !      include 'COMMON.IOUNITS'
3889 !      include 'COMMON.MD'
3890 !      include 'COMMON.CONTROL'
3891 !el local variables
3892       integer :: i
3893 !     real(kind=8) :: var,ene
3894
3895       read(inp,*) nset,nfrag,npair,nfrag_back
3896
3897 !el from module energy
3898 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
3899       if(.not.allocated(wfrag_back)) then
3900         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3901         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3902
3903         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
3904         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
3905
3906         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
3907         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
3908       endif
3909
3910       if(me.eq.king.or..not.out1file) &
3911        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
3912         " nfrag_back",nfrag_back
3913       do iset=1,nset
3914          read(inp,*) mset(iset)
3915        do i=1,nfrag
3916          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
3917            qinfrag(i,iset)
3918          if(me.eq.king.or..not.out1file) &
3919           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
3920            ifrag(2,i,iset), qinfrag(i,iset)
3921        enddo
3922        do i=1,npair
3923         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
3924           qinpair(i,iset)
3925         if(me.eq.king.or..not.out1file) &
3926          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
3927           ipair(2,i,iset), qinpair(i,iset)
3928        enddo
3929        do i=1,nfrag_back
3930         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
3931            wfrag_back(3,i,iset),&
3932            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
3933         if(me.eq.king.or..not.out1file) &
3934          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
3935          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
3936        enddo
3937       enddo
3938       return
3939       end subroutine read_fragments
3940 !-----------------------------------------------------------------------------
3941 ! shift.F       io_csa
3942 !-----------------------------------------------------------------------------
3943       subroutine csa_read
3944
3945       use csa_data
3946 !      implicit real*8 (a-h,o-z)
3947 !      include 'DIMENSIONS'
3948 !      include 'COMMON.CSA'
3949 !      include 'COMMON.BANK'
3950 !      include 'COMMON.IOUNITS'
3951 !el local variables
3952 !     integer :: ntf,ik,iw_pdb
3953 !     real(kind=8) :: var,ene
3954
3955       open(icsa_in,file=csa_in,status="old",err=100)
3956        read(icsa_in,*) nconf
3957        read(icsa_in,*) jstart,jend
3958        read(icsa_in,*) nstmax
3959        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
3960        read(icsa_in,*) nran0,nran1,irr
3961        read(icsa_in,*) nseed
3962        read(icsa_in,*) ntotal,cut1,cut2
3963        read(icsa_in,*) estop
3964        read(icsa_in,*) icmax,irestart
3965        read(icsa_in,*) ntbankm,dele,difcut
3966        read(icsa_in,*) iref,rmscut,pnccut
3967        read(icsa_in,*) ndiff
3968       close(icsa_in)
3969
3970       return
3971
3972  100  continue
3973       return
3974       end subroutine csa_read
3975 !-----------------------------------------------------------------------------
3976       subroutine initial_write
3977
3978       use csa_data
3979 !      implicit real*8 (a-h,o-z)
3980 !      include 'DIMENSIONS'
3981 !      include 'COMMON.CSA'
3982 !      include 'COMMON.BANK'
3983 !      include 'COMMON.IOUNITS'
3984 !el local variables
3985 !     integer :: ntf,ik,iw_pdb
3986 !     real(kind=8) :: var,ene
3987
3988       open(icsa_seed,file=csa_seed,status="unknown")
3989        write(icsa_seed,*) "seed"
3990       close(31)
3991 #if defined(AIX) || defined(PGI)
3992        open(icsa_history,file=csa_history,status="unknown",&
3993         position="append")
3994 #else
3995        open(icsa_history,file=csa_history,status="unknown",&
3996         access="append")
3997 #endif
3998        write(icsa_history,*) nconf
3999        write(icsa_history,*) jstart,jend
4000        write(icsa_history,*) nstmax
4001        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4002        write(icsa_history,*) nran0,nran1,irr
4003        write(icsa_history,*) nseed
4004        write(icsa_history,*) ntotal,cut1,cut2
4005        write(icsa_history,*) estop
4006        write(icsa_history,*) icmax,irestart
4007        write(icsa_history,*) ntbankm,dele,difcut
4008        write(icsa_history,*) iref,rmscut,pnccut
4009        write(icsa_history,*) ndiff
4010
4011        write(icsa_history,*)
4012       close(icsa_history)
4013
4014       open(icsa_bank1,file=csa_bank1,status="unknown")
4015        write(icsa_bank1,*) 0
4016       close(icsa_bank1)
4017
4018       return
4019       end subroutine initial_write
4020 !-----------------------------------------------------------------------------
4021       subroutine restart_write
4022
4023       use csa_data
4024 !      implicit real*8 (a-h,o-z)
4025 !      include 'DIMENSIONS'
4026 !      include 'COMMON.IOUNITS'
4027 !      include 'COMMON.CSA'
4028 !      include 'COMMON.BANK'
4029 !el local variables
4030 !     integer :: ntf,ik,iw_pdb
4031 !     real(kind=8) :: var,ene
4032
4033 #if defined(AIX) || defined(PGI)
4034        open(icsa_history,file=csa_history,position="append")
4035 #else
4036        open(icsa_history,file=csa_history,access="append")
4037 #endif
4038        write(icsa_history,*)
4039        write(icsa_history,*) "This is restart"
4040        write(icsa_history,*)
4041        write(icsa_history,*) nconf
4042        write(icsa_history,*) jstart,jend
4043        write(icsa_history,*) nstmax
4044        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4045        write(icsa_history,*) nran0,nran1,irr
4046        write(icsa_history,*) nseed
4047        write(icsa_history,*) ntotal,cut1,cut2
4048        write(icsa_history,*) estop
4049        write(icsa_history,*) icmax,irestart
4050        write(icsa_history,*) ntbankm,dele,difcut
4051        write(icsa_history,*) iref,rmscut,pnccut
4052        write(icsa_history,*) ndiff
4053        write(icsa_history,*)
4054        write(icsa_history,*) "irestart is: ", irestart
4055
4056        write(icsa_history,*)
4057       close(icsa_history)
4058
4059       return
4060       end subroutine restart_write
4061 !-----------------------------------------------------------------------------
4062 ! test.F
4063 !-----------------------------------------------------------------------------
4064       subroutine write_pdb(npdb,titelloc,ee)
4065
4066 !      implicit real*8 (a-h,o-z)
4067 !      include 'DIMENSIONS'
4068 !      include 'COMMON.IOUNITS'
4069       character(len=50) :: titelloc1
4070       character*(*) :: titelloc
4071       character(len=3) :: zahl
4072       character(len=5) :: liczba5
4073       real(kind=8) :: ee
4074       integer :: npdb   !,ilen
4075 !el      external ilen
4076 !el local variables
4077       integer :: lenpre
4078 !     real(kind=8) :: var,ene
4079
4080       titelloc1=titelloc
4081       lenpre=ilen(prefix)
4082       if (npdb.lt.1000) then
4083        call numstr(npdb,zahl)
4084        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4085       else
4086         if (npdb.lt.10000) then
4087          write(liczba5,'(i1,i4)') 0,npdb
4088         else
4089          write(liczba5,'(i5)') npdb
4090         endif
4091         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4092       endif
4093       call pdbout(ee,titelloc1,ipdb)
4094       close(ipdb)
4095       return
4096       end subroutine write_pdb
4097 !-----------------------------------------------------------------------------
4098 ! thread.F
4099 !-----------------------------------------------------------------------------
4100       subroutine write_thread_summary
4101 ! Thread the sequence through a database of known structures
4102       use control_data, only: refstr
4103 !      use geometry
4104       use energy_data, only: n_ene_comp
4105       use compare_data
4106 !      implicit real*8 (a-h,o-z)
4107 !      include 'DIMENSIONS'
4108 #ifdef MPI
4109       use MPI_data      !include 'COMMON.INFO'
4110       include 'mpif.h'
4111 #endif
4112 !      include 'COMMON.CONTROL'
4113 !      include 'COMMON.CHAIN'
4114 !      include 'COMMON.DBASE'
4115 !      include 'COMMON.INTERACT'
4116 !      include 'COMMON.VAR'
4117 !      include 'COMMON.THREAD'
4118 !      include 'COMMON.FFIELD'
4119 !      include 'COMMON.SBRIDGE'
4120 !      include 'COMMON.HEADER'
4121 !      include 'COMMON.NAMES'
4122 !      include 'COMMON.IOUNITS'
4123 !      include 'COMMON.TIME1'
4124
4125       integer,dimension(maxthread) :: ip
4126       real(kind=8),dimension(0:n_ene) :: energia
4127 !el local variables
4128       integer :: i,j,ii,jj,ipj,ik,kk,ist
4129       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4130
4131       write (iout,'(30x,a/)') &
4132        '  *********** Summary threading statistics ************'
4133       write (iout,'(a)') 'Initial energies:'
4134       write (iout,'(a4,2x,a12,14a14,3a8)') &
4135        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4136        'RMSnat','NatCONT','NNCONT','RMS'
4137 ! Energy sort patterns
4138       do i=1,nthread
4139         ip(i)=i
4140       enddo
4141       do i=1,nthread-1
4142         enet=ener(n_ene-1,ip(i))
4143         jj=i
4144         do j=i+1,nthread
4145           if (ener(n_ene-1,ip(j)).lt.enet) then
4146             jj=j
4147             enet=ener(n_ene-1,ip(j))
4148           endif
4149         enddo
4150         if (jj.ne.i) then
4151           ipj=ip(jj)
4152           ip(jj)=ip(i)
4153           ip(i)=ipj
4154         endif
4155       enddo
4156       do ik=1,nthread
4157         i=ip(ik)
4158         ii=ipatt(1,i)
4159         ist=nres_base(2,ii)+ipatt(2,i)
4160         do kk=1,n_ene_comp
4161           energia(i)=ener0(kk,i)
4162         enddo
4163         etot=ener0(n_ene_comp+1,i)
4164         rmsnat=ener0(n_ene_comp+2,i)
4165         rms=ener0(n_ene_comp+3,i)
4166         frac=ener0(n_ene_comp+4,i)
4167         frac_nn=ener0(n_ene_comp+5,i)
4168
4169         if (refstr) then
4170         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4171         i,str_nam(ii),ist+1,&
4172         (energia(print_order(kk)),kk=1,nprint_ene),&
4173         etot,rmsnat,frac,frac_nn,rms
4174         else
4175         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4176         i,str_nam(ii),ist+1,&
4177         (energia(print_order(kk)),kk=1,nprint_ene),etot
4178         endif
4179       enddo
4180       write (iout,'(//a)') 'Final energies:'
4181       write (iout,'(a4,2x,a12,17a14,3a8)') &
4182        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4183        'RMSnat','NatCONT','NNCONT','RMS'
4184       do ik=1,nthread
4185         i=ip(ik)
4186         ii=ipatt(1,i)
4187         ist=nres_base(2,ii)+ipatt(2,i)
4188         do kk=1,n_ene_comp
4189           energia(kk)=ener(kk,ik)
4190         enddo
4191         etot=ener(n_ene_comp+1,i)
4192         rmsnat=ener(n_ene_comp+2,i)
4193         rms=ener(n_ene_comp+3,i)
4194         frac=ener(n_ene_comp+4,i)
4195         frac_nn=ener(n_ene_comp+5,i)
4196         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4197         i,str_nam(ii),ist+1,&
4198         (energia(print_order(kk)),kk=1,nprint_ene),&
4199         etot,rmsnat,frac,frac_nn,rms
4200       enddo
4201       write (iout,'(/a/)') 'IEXAM array:'
4202       write (iout,'(i5)') nexcl
4203       do i=1,nexcl
4204         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4205       enddo
4206       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4207        'Max. time for threading step ',max_time_for_thread,&
4208        'Average time for threading step: ',ave_time_for_thread
4209       return
4210       end subroutine write_thread_summary
4211 !-----------------------------------------------------------------------------
4212       subroutine write_stat_thread(ithread,ipattern,ist)
4213
4214       use energy_data, only: n_ene_comp
4215       use compare_data
4216 !      implicit real*8 (a-h,o-z)
4217 !      include "DIMENSIONS"
4218 !      include "COMMON.CONTROL"
4219 !      include "COMMON.IOUNITS"
4220 !      include "COMMON.THREAD"
4221 !      include "COMMON.FFIELD"
4222 !      include "COMMON.DBASE"
4223 !      include "COMMON.NAMES"
4224       real(kind=8),dimension(0:n_ene) :: energia
4225 !el local variables
4226       integer :: ithread,ipattern,ist,i
4227       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4228
4229 #if defined(AIX) || defined(PGI)
4230       open(istat,file=statname,position='append')
4231 #else
4232       open(istat,file=statname,access='append')
4233 #endif
4234       do i=1,n_ene_comp
4235         energia(i)=ener(i,ithread)
4236       enddo
4237       etot=ener(n_ene_comp+1,ithread)
4238       rmsnat=ener(n_ene_comp+2,ithread)
4239       rms=ener(n_ene_comp+3,ithread)
4240       frac=ener(n_ene_comp+4,ithread)
4241       frac_nn=ener(n_ene_comp+5,ithread)
4242       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4243         ithread,str_nam(ipattern),ist+1,&
4244         (energia(print_order(i)),i=1,nprint_ene),&
4245         etot,rmsnat,frac,frac_nn,rms
4246       close (istat)
4247       return
4248       end subroutine write_stat_thread
4249 !-----------------------------------------------------------------------------
4250 #endif
4251 !-----------------------------------------------------------------------------
4252       end module io_config