source/unres/io_config.f90

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