source/unres/io_config.f90

   1       module io_config
   2
   3       use names
   4       use io_units
   5       use io_base
   6       use geometry_data
   7       use geometry
   8       use control_data, only:maxterm_sccor
   9       implicit none
  10 !-----------------------------------------------------------------------------
  11 ! Max. number of residue types and parameters in expressions for
  12 ! virtual-bond angle bending potentials
  13 !      integer,parameter :: maxthetyp=3
  14 !      integer,parameter :: maxthetyp1=maxthetyp+1
  15 !     ,maxtheterm=20,
  16 !     & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
  17 !     & mmaxtheterm=maxtheterm)
  18 !-----------------------------------------------------------------------------
  19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
  20 ! and the number of terms in double torsionals
  21 !      integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
  22 !      parameter (maxtor=4,maxterm=10)
  23 !-----------------------------------------------------------------------------
  24 ! Max number of torsional terms in SCCOR
  25 !el      integer,parameter :: maxterm_sccor=6
  26 !-----------------------------------------------------------------------------
  27       character(len=1),dimension(:),allocatable :: secstruc     !(maxres)
  28 !-----------------------------------------------------------------------------
  29 !
  30 !
  31 !-----------------------------------------------------------------------------
  32       contains
  33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
  34 !-----------------------------------------------------------------------------
  35 ! bank.F    io_csa
  36 !-----------------------------------------------------------------------------
  37       subroutine write_rbank(jlee,adif,nft)
  38
  39       use csa_data
  40       use geometry_data, only: nres,rad2deg
  41 !      implicit real*8 (a-h,o-z)
  42 !      include 'DIMENSIONS'
  43 !      include 'COMMON.IOUNITS'
  44 !      include 'COMMON.CSA'
  45 !      include 'COMMON.BANK'
  46 !      include 'COMMON.CHAIN'
  47 !      include 'COMMON.GEO'
  48 !el local variables
  49       integer :: nft,i,k,j,l,jlee
  50       real(kind=8) :: adif
  51
  52       open(icsa_rbank,file=csa_rbank,status="unknown")
  53       write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
  54       do k=1,nbank
  55        write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
  56        do j=1,numch
  57         do l=2,nres-1
  58          write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
  59         enddo
  60        enddo
  61       enddo
  62       close(icsa_rbank)
  63
  64   850 format (10f8.3)
  65   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
  66               i8,i10,i2,f15.5)
  67   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
  68                   ' %NC ',0pf5.2)
  69
  70       return
  71       end subroutine write_rbank
  72 !-----------------------------------------------------------------------------
  73       subroutine read_rbank(jlee,adif)
  74
  75       use csa_data
  76       use geometry_data, only: nres,deg2rad
  77       use MPI_data
  78 !      implicit real*8 (a-h,o-z)
  79 !      include 'DIMENSIONS'
  80       include 'mpif.h'
  81 !      include 'COMMON.IOUNITS'
  82 !      include 'COMMON.CSA'
  83 !      include 'COMMON.BANK'
  84 !      include 'COMMON.CHAIN'
  85 !      include 'COMMON.GEO'
  86 !      include 'COMMON.SETUP'
  87       character(len=80) :: karta
  88 !el local variables
  89       integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
  90                  ierror,ierrcode,jlee,jleer
  91       real(kind=8) :: adif
  92
  93       open(icsa_rbank,file=csa_rbank,status="old")
  94       read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
  95       print *,jleer,nbankr,nstepr,nftr,icycler,adif
  96 !       print *, 'adif from read_rbank ',adif
  97 #ifdef MPI
  98       if(nbankr.ne.nbank) then
  99         write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
 100         ' NBANK',nbank
 101         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 102       endif
 103       if(jleer.ne.jlee) then
 104         write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
 105         ' JLEE',jlee
 106         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 107       endif
 108 #endif
 109
 110       kk=0
 111       do k=1,nbankr
 112         read (icsa_rbank,'(a80)') karta
 113         write(iout,*) "READ_RBANK: kk=",kk
 114         write(iout,*) karta
 115 !        if (index(karta,"*").gt.0) then
 116 !          write (iout,*) "***** Stars in bankr ***** k=",k,
 117 !     &      " skipped"
 118 !          do j=1,numch
 119 !            do l=2,nres-1
 120 !              read (30,850) (rdummy,i=1,4)
 121 !            enddo
 122 !          enddo
 123 !        else
 124           kk=kk+1
 125           call reada(karta,"total E",rene(kk),1.0d20)
 126           call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
 127           call reada(karta,"%NC",rpncn(kk),0.0d0)
 128           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
 129             "%NC",bpncn(kk),ibank(kk)
 130 !          read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
 131           do j=1,numch
 132             do l=2,nres-1
 133               read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
 134 !              write (iout,850) (rvar(i,l,j,kk),i=1,4)
 135               do i=1,4
 136                 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
 137               enddo
 138             enddo
 139           enddo
 140 !        endif
 141       enddo
 142 !d      write (*,*) "read_rbank ******************* kk",kk,
 143 !d     &  "nbankr",nbankr
 144       if (kk.lt.nbankr) nbankr=kk
 145 !d      do kk=1,nbankr
 146 !d          print *,"kk=",kk
 147 !d          do j=1,numch
 148 !d            do l=2,nres-1
 149 !d              write (*,850) (rvar(i,l,j,kk),i=1,4)
 150 !d            enddo
 151 !d          enddo
 152 !d      enddo
 153       close(icsa_rbank)
 154
 155   850 format (10f8.3)
 156   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 157   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
 158
 159       return
 160       end subroutine read_rbank
 161 !-----------------------------------------------------------------------------
 162       subroutine write_bank(jlee,nft)
 163
 164       use csa_data
 165       use control_data, only: vdisulf
 166       use geometry_data, only: nres,rad2deg
 167 !      implicit real*8 (a-h,o-z)
 168 !      include 'DIMENSIONS'
 169 !      include 'COMMON.IOUNITS'
 170 !      include 'COMMON.CSA'
 171 !      include 'COMMON.BANK'
 172 !      include 'COMMON.CHAIN'
 173 !      include 'COMMON.GEO'
 174 !      include 'COMMON.SBRIDGE'
 175 !     include 'COMMON.CONTROL'
 176       character(len=7) :: chtmp
 177       character(len=40) :: chfrm
 178 !el      external ilen
 179 !el local variables
 180       integer :: nft,k,l,i,j,jlee
 181
 182       open(icsa_bank,file=csa_bank,status="unknown")
 183       write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
 184       write (icsa_bank,902) nglob_csa, eglob_csa
 185       open (igeom,file=intname,status='UNKNOWN')
 186       do k=1,nbank
 187        write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 188        if (vdisulf) write (icsa_bank,'(101i4)') &
 189           bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
 190        do j=1,numch
 191         do l=2,nres-1
 192          write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 193         enddo
 194        enddo
 195        if (bvar_nss(k).le.9) then
 196          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 197            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 198        else
 199          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 200            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 201          write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
 202                                       bvar_ss(2,i,k),i=10,bvar_nss(k))
 203        endif
 204        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 205        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 206        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 207        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 208       enddo
 209       close(icsa_bank)
 210       close(igeom)
 211
 212       if (nstep/200.gt.ilastnstep) then
 213
 214        ilastnstep=(ilastnstep+1)*1.5
 215        write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
 216        write(chtmp,chfrm) nstep
 217        open(icsa_int,file=prefix(:ilen(prefix)) &
 218                //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
 219        do k=1,nbank
 220         if (bvar_nss(k).le.9) then
 221          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 222         bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 223         else
 224          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 225            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 226          write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
 227                                 bvar_ss(2,i,k),i=10,bvar_nss(k))
 228         endif
 229         write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 230         write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 231         write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 232         write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 233        enddo
 234        close(icsa_int)
 235       endif
 236
 237
 238   200 format (8f10.4)
 239   850 format (10f8.3)
 240   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
 241               i8,i10,i2,f15.5)
 242   902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
 243   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 244               ' %NC ',0pf5.2,i5)
 245
 246       return
 247       end subroutine write_bank
 248 !-----------------------------------------------------------------------------
 249       subroutine write_bank_reminimized(jlee,nft)
 250
 251       use csa_data
 252       use geometry_data, only: nres,rad2deg
 253       use energy_data
 254 !      implicit real*8 (a-h,o-z)
 255 !      include 'DIMENSIONS'
 256 !      include 'COMMON.IOUNITS'
 257 !      include 'COMMON.CSA'
 258 !      include 'COMMON.BANK'
 259 !      include 'COMMON.CHAIN'
 260 !      include 'COMMON.GEO'
 261 !      include 'COMMON.SBRIDGE'
 262 !el local variables
 263       integer :: nft,i,l,j,k,jlee
 264
 265       open(icsa_bank_reminimized,file=csa_bank_reminimized,&
 266         status="unknown")
 267       write (icsa_bank_reminimized,900) &
 268         jlee,nbank,nstep,nft,icycle,cutdif
 269       open (igeom,file=intname,status='UNKNOWN')
 270       do k=1,nbank
 271        write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
 272         bpncn(k),ibank(k)
 273        do j=1,numch
 274         do l=2,nres-1
 275          write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 276         enddo
 277        enddo
 278        if (nss.le.9) then
 279          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 280            nss,(ihpb(i),jhpb(i),i=1,nss)
 281        else
 282          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 283            nss,(ihpb(i),jhpb(i),i=1,9)
 284          write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
 285        endif
 286        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 287        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 288        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 289        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 290       enddo
 291       close(icsa_bank_reminimized)
 292       close(igeom)
 293
 294   200 format (8f10.4)
 295   850 format (10f8.3)
 296   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
 297               i8,i10,i2,f15.5)
 298   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 299                ' %NC ',0pf5.2,i5)
 300
 301       return
 302       end subroutine write_bank_reminimized
 303 !-----------------------------------------------------------------------------
 304       subroutine read_bank(jlee,nft,cutdifr)
 305
 306       use csa_data
 307       use control_data, only: vdisulf
 308       use geometry_data, only: nres,deg2rad
 309       use energy_data
 310 !      implicit real*8 (a-h,o-z)
 311 !      include 'DIMENSIONS'
 312 !      include 'COMMON.IOUNITS'
 313 !      include 'COMMON.CSA'
 314 !      include 'COMMON.BANK'
 315 !      include 'COMMON.CHAIN'
 316 !      include 'COMMON.GEO'
 317 !      include 'COMMON.CONTROL'
 318 !      include 'COMMON.SBRIDGE'
 319       character(len=80) :: karta
 320 !      integer ilen
 321 !el      external ilen
 322 !el local variables
 323       integer :: nft,kk,k,l,i,j,jlee
 324       real(kind=8) :: cutdifr
 325
 326       open(icsa_bank,file=csa_bank,status="old")
 327        read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
 328        read (icsa_bank,902) nglob_csa, eglob_csa
 329 !      if(jleer.ne.jlee) then
 330 !        write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
 331 !    &   ' JLEE',jlee
 332 !        call mpi_abort(mpi_comm_world,ierror,ierrcode)
 333 !      endif
 334
 335       kk=0
 336       do k=1,nbank
 337         read (icsa_bank,'(a80)') karta
 338         write(iout,*) "READ_BANK: kk=",kk
 339         write(iout,*) karta
 340 !        if (index(karta,"*").gt.0) then
 341 !          write (iout,*) "***** Stars in bank ***** k=",k,
 342 !     &      " skipped"
 343 !          do j=1,numch
 344 !            do l=2,nres-1
 345 !              read (33,850) (rdummy,i=1,4)
 346 !            enddo
 347 !          enddo
 348 !        else
 349           kk=kk+1
 350           call reada(karta,"total E",bene(kk),1.0d20)
 351           call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
 352           call reada(karta,"%NC",bpncn(kk),0.0d0)
 353           read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
 354           goto 112
 355   111     ibank(kk)=0
 356   112     continue
 357           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
 358             "%NC",bpncn(kk),ibank(kk)
 359 !          read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 360           if (vdisulf) then
 361             read (icsa_bank,'(101i4)') &
 362               bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 363             bvar_ns(kk)=ns-2*bvar_nss(kk)
 364             write(iout,*) 'read SSBOND',bvar_nss(kk),&
 365                           ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 366 !d          write(iout,*) 'read CYS #free ', bvar_ns(kk)
 367             l=0
 368             do i=1,ns
 369              j=1
 370              do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
 371                        iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
 372                        j.le.bvar_nss(kk))
 373                 j=j+1
 374              enddo
 375              if (j.gt.bvar_nss(kk)) then
 376                l=l+1
 377                bvar_s(l,kk)=iss(i)
 378              endif
 379             enddo
 380 !d            write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
 381           endif
 382           do j=1,numch
 383             do l=2,nres-1
 384               read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
 385 !              write (iout,850) (bvar(i,l,j,kk),i=1,4)
 386               do i=1,4
 387                 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
 388               enddo ! l
 389             enddo ! l
 390           enddo ! j
 391 !        endif
 392       enddo ! k
 393
 394       if (kk.lt.nbank) nbank=kk
 395 !d      write (*,*) "read_bank ******************* kk",kk,
 396 !d     &  "nbank",nbank
 397 !d      do kk=1,nbank
 398 !d          print *,"kk=",kk
 399 !d          do j=1,numch
 400 !d            do l=2,nres-1
 401 !d              write (*,850) (bvar(i,l,j,kk),i=1,4)
 402 !d            enddo
 403 !d          enddo
 404 !d      enddo
 405
 406 !       do k=1,nbank
 407 !        read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 408 !        do j=1,numch
 409 !         do l=2,nres-1
 410 !          read (33,850) (bvar(i,l,j,k),i=1,4)
 411 !          do i=1,4
 412 !           bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
 413 !          enddo
 414 !         enddo
 415 !        enddo
 416 !       enddo
 417       close(icsa_bank)
 418
 419   850 format (10f8.3)
 420   952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
 421   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 422   902 format (1x,11x,i4,12x,1pe14.5)
 423   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
 424
 425       return
 426       end subroutine read_bank
 427 !-----------------------------------------------------------------------------
 428       subroutine write_bank1(jlee)
 429
 430       use csa_data
 431       use geometry_data, only: nres,rad2deg
 432 !      implicit real*8 (a-h,o-z)
 433 !      include 'DIMENSIONS'
 434 !      include 'COMMON.IOUNITS'
 435 !      include 'COMMON.CSA'
 436 !      include 'COMMON.BANK'
 437 !      include 'COMMON.CHAIN'
 438 !      include 'COMMON.GEO'
 439 !el local variables
 440       integer :: k,i,l,j,jlee
 441
 442 #if defined(AIX) || defined(PGI)
 443       open(icsa_bank1,file=csa_bank1,position="append")
 444 #else
 445       open(icsa_bank1,file=csa_bank1,access="append")
 446 #endif
 447       write (icsa_bank1,900) jlee,nbank,nstep,cutdif
 448       do k=1,nbank
 449        write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 450        do j=1,numch
 451         do l=2,nres-1
 452          write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 453         enddo
 454        enddo
 455       enddo
 456       close(icsa_bank1)
 457   850 format (10f8.3)
 458   900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
 459   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 460               ' %NC ',0pf5.2,i5)
 461
 462       return
 463       end subroutine write_bank1
 464 !-----------------------------------------------------------------------------
 465 ! cartprint.f
 466 !-----------------------------------------------------------------------------
 467 !      subroutine cartprint
 468
 469 !      use geometry_data, only: c
 470 !      use energy_data, only: itype
 471 !      implicit real*8 (a-h,o-z)
 472 !      include 'DIMENSIONS'
 473 !      include 'COMMON.CHAIN'
 474 !      include 'COMMON.INTERACT'
 475 !      include 'COMMON.NAMES'
 476 !      include 'COMMON.IOUNITS'
 477 !      integer :: i
 478
 479 !     write (iout,100)
 480 !      do i=1,nres
 481 !        write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
 482 !          c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
 483 !      enddo
 484 !  100 format (//'              alpha-carbon coordinates       ',&
 485 !                '     centroid coordinates'/ &
 486 !                '       ', 6X,'X',11X,'Y',11X,'Z',&
 487 !                                10X,'X',11X,'Y',11X,'Z')
 488 !  110 format (a,'(',i3,')',6f12.5)
 489 !      return
 490 !      end subroutine cartprint
 491 !-----------------------------------------------------------------------------
 492 ! dihed_cons.F
 493 !-----------------------------------------------------------------------------
 494       subroutine secstrp2dihc
 495
 496       use geometry_data
 497       use energy_data
 498 !      implicit real*8 (a-h,o-z)
 499 !      include 'DIMENSIONS'
 500 !      include 'COMMON.GEO'
 501 !      include 'COMMON.BOUNDS'
 502 !      include 'COMMON.CHAIN'
 503 !      include 'COMMON.TORCNSTR'
 504 !      include 'COMMON.IOUNITS'
 505 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 506 !el      COMMON/SECONDARYS/secstruc
 507       character(len=80) :: line
 508       logical :: errflag
 509 !el      external ilen
 510
 511 !el  local variables
 512       integer :: i,ii,lenpre
 513
 514       allocate(secstruc(nres))
 515
 516 !dr      call getenv_loc('SECPREDFIL',secpred)
 517       lenpre=ilen(prefix)
 518       secpred=prefix(:lenpre)//'.spred'
 519
 520 #if defined(WINIFL) || defined(WINPGI)
 521       open(isecpred,file=secpred,status='old',readonly,shared)
 522 #elif (defined CRAY) || (defined AIX)
 523       open(isecpred,file=secpred,status='old',action='read')
 524 #elif (defined G77)
 525       open(isecpred,file=secpred,status='old')
 526 #else
 527       open(isecpred,file=secpred,status='old',action='read')
 528 #endif
 529 ! read secondary structure prediction from JPRED here!
 530 !      read(isecpred,'(A80)',err=100,end=100) line
 531 !      read(line,'(f10.3)',err=110) ftors
 532        read(isecpred,'(f10.3)',err=110) ftors
 533
 534       write (iout,*) 'FTORS factor =',ftors
 535 ! initialize secstruc to any
 536        do i=1,nres
 537         secstruc(i) ='-'
 538       enddo
 539       ndih_constr=0
 540       ndih_nconstr=0
 541
 542       call read_secstr_pred(isecpred,iout,errflag)
 543       if (errflag) then
 544          write(iout,*)'There is a problem with the list of secondary-',&
 545            'structure prediction'
 546          goto 100
 547       endif
 548 ! 8/13/98 Set limits to generating the dihedral angles
 549       do i=1,nres
 550         phibound(1,i)=-pi
 551         phibound(2,i)=pi
 552       enddo
 553
 554       ii=0
 555       do i=1,nres
 556         if ( secstruc(i) .eq. 'H') then
 557 ! Helix restraints for this residue
 558            ii=ii+1
 559            idih_constr(ii)=i
 560            phi0(ii) = 45.0D0*deg2rad
 561            drange(ii)= 5.0D0*deg2rad
 562            phibound(1,i) = phi0(ii)-drange(ii)
 563            phibound(2,i) = phi0(ii)+drange(ii)
 564         else if (secstruc(i) .eq. 'E') then
 565 ! strand restraints for this residue
 566            ii=ii+1
 567            idih_constr(ii)=i
 568            phi0(ii) = 180.0D0*deg2rad
 569            drange(ii)= 5.0D0*deg2rad
 570            phibound(1,i) = phi0(ii)-drange(ii)
 571            phibound(2,i) = phi0(ii)+drange(ii)
 572         else
 573 ! no restraints for this residue
 574            ndih_nconstr=ndih_nconstr+1
 575            idih_nconstr(ndih_nconstr)=i
 576         endif
 577       enddo
 578       ndih_constr=ii
 579 !      deallocate(secstruc)
 580       return
 581 100   continue
 582       write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
 583 !      deallocate(secstruc)
 584       return
 585 110   continue
 586       write(iout,'(A20)')'Error reading FTORS'
 587 !      deallocate(secstruc)
 588       return
 589       end subroutine secstrp2dihc
 590 !-----------------------------------------------------------------------------
 591       subroutine read_secstr_pred(jin,jout,errors)
 592
 593 !      implicit real*8 (a-h,o-z)
 594 !      INCLUDE 'DIMENSIONS'
 595 !      include 'COMMON.IOUNITS'
 596 !      include 'COMMON.CHAIN'
 597 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 598 !el      COMMON/SECONDARYS/secstruc
 599 !el      EXTERNAL ILEN
 600       character(len=80) :: line,line1   !,ucase
 601       logical :: errflag,errors,blankline
 602
 603 !el  local variables
 604       integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
 605             length_of_chain
 606       errors=.false.
 607       read (jin,'(a)') line
 608       write (jout,'(2a)') '> ',line(1:78)
 609       line1=ucase(line)
 610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
 611 ! correspond to the end-groups.  ADD to the secondary structure prediction "-" for the
 612 ! end-groups in the input file "*.spred"
 613
 614       iseq=1
 615       do while (index(line1,'$END').eq.0)
 616 ! Override commented lines.
 617          ipos=1
 618          blankline=.false.
 619          do while (.not.blankline)
 620             line1=' '
 621             call mykey(line,line1,ipos,blankline,errflag)
 622             if (errflag) write (jout,'(2a)') &
 623        'Error when reading sequence in line: ',line
 624             errors=errors .or. errflag
 625             if (.not. blankline .and. .not. errflag) then
 626                ipos1=2
 627                iend=ilen(line1)
 628 !el               if (iseq.le.maxres) then
 629                   if (line1(1:1).eq.'-' ) then
 630                      secstruc(iseq)=line1(1:1)
 631                   else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
 632                             ( ucase(line1(1:1)).eq.'H' ) ) then
 633                      secstruc(iseq)=ucase(line1(1:1))
 634                   else
 635                      errors=.true.
 636                      write (jout,1010) line1(1:1), iseq
 637                      goto 80
 638                   endif
 639 !el               else
 640 !el                  errors=.true.
 641 !el                  write (jout,1000) iseq,maxres
 642 !el                  goto 80
 643 !el               endif
 644                do while (ipos1.le.iend)
 645
 646                   iseq=iseq+1
 647                   il=1
 648                   ipos1=ipos1+1
 649 !el                  if (iseq.le.maxres) then
 650                      if (line1(ipos1-1:ipos1-1).eq.'-' ) then
 651                         secstruc(iseq)=line1(ipos1-1:ipos1-1)
 652                      else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
 653                            (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
 654                         secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
 655                      else
 656                         errors=.true.
 657                         write (jout,1010) line1(ipos1-1:ipos1-1), iseq
 658                         goto 80
 659                      endif
 660 !el                  else
 661 !el                     errors=.true.
 662 !el                     write (jout,1000) iseq,maxres
 663 !el                     goto 80
 664 !el                  endif
 665                enddo
 666                iseq=iseq+1
 667             endif
 668          enddo
 669          read (jin,'(a)') line
 670          write (jout,'(2a)') '> ',line(1:78)
 671          line1=ucase(line)
 672       enddo
 673
 674 !d    write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
 675
 676 !d check whether the found length of the chain is correct.
 677       length_of_chain=iseq-1
 678       if (length_of_chain .ne. nres) then
 679 !        errors=.true.
 680         write (jout,'(a,i4,a,i4,a)') &
 681        'Error: the number of labels specified in $SEC_STRUC_PRED (' &
 682        ,length_of_chain,') does not match with the number of residues (' &
 683        ,nres,').'
 684       endif
 685    80 continue
 686
 687  1000 format('Error - the number of residues (',i4,&
 688        ') has exceeded maximum (',i4,').')
 689  1010 format ('Error - unrecognized secondary structure label',a4,&
 690        ' in position',i4)
 691       return
 692       end subroutine read_secstr_pred
 693 !#endif
 694 !-----------------------------------------------------------------------------
 695 ! parmread.F
 696 !-----------------------------------------------------------------------------
 697       subroutine parmread
 698
 699       use geometry_data
 700       use energy_data
 701       use control_data, only:maxterm !,maxtor
 702       use MD_data
 703       use MPI_data
 704 !el      use map_data
 705       use control, only: getenv_loc
 706 !
 707 ! Read the parameters of the probability distributions of the virtual-bond
 708 ! valence angles and the side chains and energy parameters.
 709 !
 710 ! Important! Energy-term weights ARE NOT read here; they are read from the
 711 ! main input file instead, because NO defaults have yet been set for these
 712 ! parameters.
 713 !
 714 !      implicit real*8 (a-h,o-z)
 715 !      include 'DIMENSIONS'
 716 #ifdef MPI
 717       include "mpif.h"
 718       integer :: IERROR
 719 #endif
 720 !      include 'COMMON.IOUNITS'
 721 !      include 'COMMON.CHAIN'
 722 !      include 'COMMON.INTERACT'
 723 !      include 'COMMON.GEO'
 724 !      include 'COMMON.LOCAL'
 725 !      include 'COMMON.TORSION'
 726 !      include 'COMMON.SCCOR'
 727 !      include 'COMMON.SCROT'
 728 !      include 'COMMON.FFIELD'
 729 !      include 'COMMON.NAMES'
 730 !      include 'COMMON.SBRIDGE'
 731 !      include 'COMMON.MD'
 732 !      include 'COMMON.SETUP'
 733       character(len=1) :: t1,t2,t3
 734       character(len=1) :: onelett(4) = (/"G","A","P","D"/)
 735       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
 736       logical :: lprint,LaTeX
 737       real(kind=8),dimension(3,3,maxlob) :: blower      !(3,3,maxlob)
 738       real(kind=8),dimension(13) :: b
 739       character(len=3) :: lancuch       !,ucase
 740 !el  local variables
 741       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
 742       integer :: maxinter,junk,kk,ii
 743       real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
 744                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
 745                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
 746                 res1,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),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),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),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),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),&
 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),&
 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),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),&
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),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),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),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),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),restyp(j),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)=ntyp1
2389           itype(ires_old-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)=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)=10
2460             else
2461               itype(ires)=rescode(ires,res,0)
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),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)
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)
2501 !        if (itype(i).eq.ntyp1) then
2502 !          write (iout,*) "dummy",i,itype(i)
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).eq.ntyp1) then
2510          if (itype(i+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).eq.ntyp1)=true
2513 ! second dummy atom is conected to next chain itype(i+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).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).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).ne.10) then
2573         nres=nres+1
2574         itype(nres)=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).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)),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),restyp(itype(ires)),(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).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)),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       write (iout,*) "with_theta_constr ",with_theta_constr
2975       AFMlog=(index(controlcard,'AFM'))
2976       selfguide=(index(controlcard,'SELFGUIDE'))
2977       print *,'AFMlog',AFMlog,selfguide,"KUPA"
2978       call readi(controlcard,'GENCONSTR',genconstr,0)
2979       call reada(controlcard,'BOXX',boxxsize,100.0d0)
2980       call reada(controlcard,'BOXY',boxysize,100.0d0)
2981       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2982       call readi(controlcard,'TUBEMOD',tubemode,0)
2983       write (iout,*) TUBEmode,"TUBEMODE"
2984       if (TUBEmode.gt.0) then
2985        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
2986        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
2987        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
2988        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
2989        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
2990        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
2991        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
2992        buftubebot=bordtubebot+tubebufthick
2993        buftubetop=bordtubetop-tubebufthick
2994       endif
2995
2996 ! CUTOFFF ON ELECTROSTATICS
2997       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2998       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2999       write(iout,*) "R_CUT_ELE=",r_cut_ele
3000 ! Lipidic parameters
3001       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3002       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3003       if (lipthick.gt.0.0d0) then
3004        bordliptop=(boxzsize+lipthick)/2.0
3005        bordlipbot=bordliptop-lipthick
3006       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3007       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3008       buflipbot=bordlipbot+lipbufthick
3009       bufliptop=bordliptop-lipbufthick
3010       if ((lipbufthick*2.0d0).gt.lipthick) &
3011        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3012       endif !lipthick.gt.0
3013       write(iout,*) "bordliptop=",bordliptop
3014       write(iout,*) "bordlipbot=",bordlipbot
3015       write(iout,*) "bufliptop=",bufliptop
3016       write(iout,*) "buflipbot=",buflipbot
3017       write (iout,*) "SHIELD MODE",shield_mode
3018
3019 !C-------------------------
3020       minim=(index(controlcard,'MINIMIZE').gt.0)
3021       dccart=(index(controlcard,'CART').gt.0)
3022       overlapsc=(index(controlcard,'OVERLAP').gt.0)
3023       overlapsc=.not.overlapsc
3024       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3025       searchsc=.not.searchsc
3026       sideadd=(index(controlcard,'SIDEADD').gt.0)
3027       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3028       outpdb=(index(controlcard,'PDBOUT').gt.0)
3029       outmol2=(index(controlcard,'MOL2OUT').gt.0)
3030       pdbref=(index(controlcard,'PDBREF').gt.0)
3031       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3032       indpdb=index(controlcard,'PDBSTART')
3033       extconf=(index(controlcard,'EXTCONF').gt.0)
3034       call readi(controlcard,'IPRINT',iprint,0)
3035       call readi(controlcard,'MAXGEN',maxgen,10000)
3036       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3037       call readi(controlcard,"KDIAG",kdiag,0)
3038       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3039       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3040        write (iout,*) "RESCALE_MODE",rescale_mode
3041       split_ene=index(controlcard,'SPLIT_ENE').gt.0
3042       if (index(controlcard,'REGULAR').gt.0.0D0) then
3043         call reada(controlcard,'WEIDIS',weidis,0.1D0)
3044         modecalc=1
3045         refstr=.true.
3046       endif
3047       if (index(controlcard,'CHECKGRAD').gt.0) then
3048         modecalc=5
3049         if (index(controlcard,'CART').gt.0) then
3050           icheckgrad=1
3051         elseif (index(controlcard,'CARINT').gt.0) then
3052           icheckgrad=2
3053         else
3054           icheckgrad=3
3055         endif
3056       elseif (index(controlcard,'THREAD').gt.0) then
3057         modecalc=2
3058         call readi(controlcard,'THREAD',nthread,0)
3059         if (nthread.gt.0) then
3060           call reada(controlcard,'WEIDIS',weidis,0.1D0)
3061         else
3062           if (fg_rank.eq.0) &
3063           write (iout,'(a)')'A number has to follow the THREAD keyword.'
3064           stop 'Error termination in Read_Control.'
3065         endif
3066       else if (index(controlcard,'MCMA').gt.0) then
3067         modecalc=3
3068       else if (index(controlcard,'MCEE').gt.0) then
3069         modecalc=6
3070       else if (index(controlcard,'MULTCONF').gt.0) then
3071         modecalc=4
3072       else if (index(controlcard,'MAP').gt.0) then
3073         modecalc=7
3074         call readi(controlcard,'MAP',nmap,0)
3075       else if (index(controlcard,'CSA').gt.0) then
3076         modecalc=8
3077 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3078 !rc
3079 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3080 !rc
3081 !rc        modecalc=9
3082 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3083 !fmc        modecalc=10
3084       else if (index(controlcard,'SOFTREG').gt.0) then
3085         modecalc=11
3086       else if (index(controlcard,'CHECK_BOND').gt.0) then
3087         modecalc=-1
3088       else if (index(controlcard,'TEST').gt.0) then
3089         modecalc=-2
3090       else if (index(controlcard,'MD').gt.0) then
3091         modecalc=12
3092       else if (index(controlcard,'RE ').gt.0) then
3093         modecalc=14
3094       endif
3095
3096       lmuca=index(controlcard,'MUCA').gt.0
3097       call readi(controlcard,'MUCADYN',mucadyn,0)
3098       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3099       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3100        then
3101        write (iout,*) 'MUCADYN=',mucadyn
3102        write (iout,*) 'MUCASMOOTH=',muca_smooth
3103       endif
3104
3105       iscode=index(controlcard,'ONE_LETTER')
3106       indphi=index(controlcard,'PHI')
3107       indback=index(controlcard,'BACK')
3108       iranconf=index(controlcard,'RAND_CONF')
3109       i2ndstr=index(controlcard,'USE_SEC_PRED')
3110       gradout=index(controlcard,'GRADOUT').gt.0
3111       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3112       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3113       if (me.eq.king .or. .not.out1file ) &
3114         write (iout,*) "DISTCHAINMAX",distchainmax
3115
3116       if(me.eq.king.or..not.out1file) &
3117        write (iout,'(2a)') diagmeth(kdiag),&
3118         ' routine used to diagonalize matrices.'
3119       if (shield_mode.gt.0) then
3120       pi=3.141592d0
3121 !C VSolvSphere the volume of solving sphere
3122 !C      print *,pi,"pi"
3123 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3124 !C there will be no distinction between proline peptide group and normal peptide
3125 !C group in case of shielding parameters
3126       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3127       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3128       write (iout,*) VSolvSphere,VSolvSphere_div
3129 !C long axis of side chain
3130       do i=1,ntyp
3131       long_r_sidechain(i)=vbldsc0(1,i)
3132       short_r_sidechain(i)=sigma0(i)
3133       write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3134          sigma0(i)
3135       enddo
3136       buff_shield=1.0d0
3137       endif
3138       return
3139       end subroutine read_control
3140 !-----------------------------------------------------------------------------
3141       subroutine read_REMDpar
3142 !
3143 ! Read REMD settings
3144 !
3145 !       use control
3146 !       use energy
3147 !       use geometry
3148       use REMD_data
3149       use MPI_data
3150       use control_data, only:out1file
3151 !      implicit real*8 (a-h,o-z)
3152 !      include 'DIMENSIONS'
3153 !      include 'COMMON.IOUNITS'
3154 !      include 'COMMON.TIME1'
3155 !      include 'COMMON.MD'
3156       use MD_data
3157 !el #ifndef LANG0
3158 !el      include 'COMMON.LANGEVIN'
3159 !el #else
3160 !el      include 'COMMON.LANGEVIN.lang0'
3161 !el #endif
3162 !      include 'COMMON.INTERACT'
3163 !      include 'COMMON.NAMES'
3164 !      include 'COMMON.GEO'
3165 !      include 'COMMON.REMD'
3166 !      include 'COMMON.CONTROL'
3167 !      include 'COMMON.SETUP'
3168 !      character(len=80) :: ucase
3169       character(len=320) :: controlcard
3170       character(len=3200) :: controlcard1
3171       integer :: iremd_m_total
3172 !el local variables
3173       integer :: i
3174 !     real(kind=8) :: var,ene
3175
3176       if(me.eq.king.or..not.out1file) &
3177        write (iout,*) "REMD setup"
3178
3179       call card_concat(controlcard,.true.)
3180       call readi(controlcard,"NREP",nrep,3)
3181       call readi(controlcard,"NSTEX",nstex,1000)
3182       call reada(controlcard,"RETMIN",retmin,10.0d0)
3183       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3184       mremdsync=(index(controlcard,'SYNC').gt.0)
3185       call readi(controlcard,"NSYN",i_sync_step,100)
3186       restart1file=(index(controlcard,'REST1FILE').gt.0)
3187       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3188       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3189       if(max_cache_traj_use.gt.max_cache_traj) &
3190                 max_cache_traj_use=max_cache_traj
3191       if(me.eq.king.or..not.out1file) then
3192 !d       if (traj1file) then
3193 !rc caching is in testing - NTWX is not ignored
3194 !d        write (iout,*) "NTWX value is ignored"
3195 !d        write (iout,*) "  trajectory is stored to one file by master"
3196 !d        write (iout,*) "  before exchange at NSTEX intervals"
3197 !d       endif
3198        write (iout,*) "NREP= ",nrep
3199        write (iout,*) "NSTEX= ",nstex
3200        write (iout,*) "SYNC= ",mremdsync
3201        write (iout,*) "NSYN= ",i_sync_step
3202        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3203       endif
3204       remd_tlist=.false.
3205       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3206       if (index(controlcard,'TLIST').gt.0) then
3207          remd_tlist=.true.
3208          call card_concat(controlcard1,.true.)
3209          read(controlcard1,*) (remd_t(i),i=1,nrep)
3210          if(me.eq.king.or..not.out1file) &
3211           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3212       endif
3213       remd_mlist=.false.
3214       if (index(controlcard,'MLIST').gt.0) then
3215          remd_mlist=.true.
3216          call card_concat(controlcard1,.true.)
3217          read(controlcard1,*) (remd_m(i),i=1,nrep)
3218          if(me.eq.king.or..not.out1file) then
3219           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3220           iremd_m_total=0
3221           do i=1,nrep
3222            iremd_m_total=iremd_m_total+remd_m(i)
3223           enddo
3224           write (iout,*) 'Total number of replicas ',iremd_m_total
3225          endif
3226       endif
3227       if(me.eq.king.or..not.out1file) &
3228        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3229       return
3230       end subroutine read_REMDpar
3231 !-----------------------------------------------------------------------------
3232       subroutine read_MDpar
3233 !
3234 ! Read MD settings
3235 !
3236       use control_data, only: r_cut,rlamb,out1file
3237       use energy_data
3238       use geometry_data, only: pi
3239       use MPI_data
3240 !      implicit real*8 (a-h,o-z)
3241 !      include 'DIMENSIONS'
3242 !      include 'COMMON.IOUNITS'
3243 !      include 'COMMON.TIME1'
3244 !      include 'COMMON.MD'
3245       use MD_data
3246 !el #ifndef LANG0
3247 !el      include 'COMMON.LANGEVIN'
3248 !el #else
3249 !el      include 'COMMON.LANGEVIN.lang0'
3250 !el #endif
3251 !      include 'COMMON.INTERACT'
3252 !      include 'COMMON.NAMES'
3253 !      include 'COMMON.GEO'
3254 !      include 'COMMON.SETUP'
3255 !      include 'COMMON.CONTROL'
3256 !      include 'COMMON.SPLITELE'
3257 !      character(len=80) :: ucase
3258       character(len=320) :: controlcard
3259 !el local variables
3260       integer :: i
3261       real(kind=8) :: eta
3262
3263       call card_concat(controlcard,.true.)
3264       call readi(controlcard,"NSTEP",n_timestep,1000000)
3265       call readi(controlcard,"NTWE",ntwe,100)
3266       call readi(controlcard,"NTWX",ntwx,1000)
3267       call reada(controlcard,"DT",d_time,1.0d-1)
3268       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3269       call reada(controlcard,"DAMAX",damax,1.0d1)
3270       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3271       call readi(controlcard,"LANG",lang,0)
3272       RESPA = index(controlcard,"RESPA") .gt. 0
3273       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3274       ntime_split0=ntime_split
3275       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3276       ntime_split0=ntime_split
3277       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3278       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3279       rest = index(controlcard,"REST").gt.0
3280       tbf = index(controlcard,"TBF").gt.0
3281       usampl = index(controlcard,"USAMPL").gt.0
3282       mdpdb = index(controlcard,"MDPDB").gt.0
3283       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3284       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3285       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3286       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3287       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3288       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3289       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3290       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3291       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3292       large = index(controlcard,"LARGE").gt.0
3293       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3294       rattle = index(controlcard,"RATTLE").gt.0
3295 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3296       nset=0
3297       if(usampl) then
3298         call read_fragments
3299       endif
3300
3301       if(me.eq.king.or..not.out1file) then
3302        write (iout,*)
3303        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3304        write (iout,*)
3305        write (iout,'(a)') "The units are:"
3306        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3307        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3308         " acceleration: angstrom/(48.9 fs)**2"
3309        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3310        write (iout,*)
3311        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3312        write (iout,'(a60,f10.5,a)') &
3313         "Initial time step of numerical integration:",d_time,&
3314         " natural units"
3315        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3316        if (RESPA) then
3317         write (iout,'(2a,i4,a)') &
3318           "A-MTS algorithm used; initial time step for fast-varying",&
3319           " short-range forces split into",ntime_split," steps."
3320         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3321          r_cut," lambda",rlamb
3322        endif
3323        write (iout,'(2a,f10.5)') &
3324         "Maximum acceleration threshold to reduce the time step",&
3325         "/increase split number:",damax
3326        write (iout,'(2a,f10.5)') &
3327         "Maximum predicted energy drift to reduce the timestep",&
3328         "/increase split number:",edriftmax
3329        write (iout,'(a60,f10.5)') &
3330        "Maximum velocity threshold to reduce velocities:",dvmax
3331        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3332        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3333        if (rattle) write (iout,'(a60)') &
3334         "Rattle algorithm used to constrain the virtual bonds"
3335       endif
3336       reset_fricmat=1000
3337       if (lang.gt.0) then
3338         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3339         call reada(controlcard,"RWAT",rwat,1.4d0)
3340         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3341         surfarea=index(controlcard,"SURFAREA").gt.0
3342         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3343         if(me.eq.king.or..not.out1file)then
3344          write (iout,'(/a,$)') "Langevin dynamics calculation"
3345          if (lang.eq.1) then
3346           write (iout,'(a/)') &
3347             " with direct integration of Langevin equations"
3348          else if (lang.eq.2) then
3349           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3350          else if (lang.eq.3) then
3351           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3352          else if (lang.eq.4) then
3353           write (iout,'(a/)') " in overdamped mode"
3354          else
3355           write (iout,'(//a,i5)') &
3356             "=========== ERROR: Unknown Langevin dynamics mode:",lang
3357           stop
3358          endif
3359          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3360          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3361          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3362          write (iout,'(a60,f10.5)') &
3363          "Scaling factor of the friction forces:",scal_fric
3364          if (surfarea) write (iout,'(2a,i10,a)') &
3365            "Friction coefficients will be scaled by solvent-accessible",&
3366            " surface area every",reset_fricmat," steps."
3367         endif
3368 ! Calculate friction coefficients and bounds of stochastic forces
3369         eta=6*pi*cPoise*etawat
3370         if(me.eq.king.or..not.out1file) &
3371          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3372           eta
3373         gamp=scal_fric*(pstok+rwat)*eta
3374         stdfp=dsqrt(2*Rb*t_bath/d_time)
3375         allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3376         do i=1,ntyp
3377           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3378           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3379         enddo
3380         if(me.eq.king.or..not.out1file)then
3381          write (iout,'(/2a/)') &
3382          "Radii of site types and friction coefficients and std's of",&
3383          " stochastic forces of fully exposed sites"
3384          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3385          do i=1,ntyp
3386           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3387            gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3388          enddo
3389         endif
3390       else if (tbf) then
3391         if(me.eq.king.or..not.out1file)then
3392          write (iout,'(a)') "Berendsen bath calculation"
3393          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3394          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3395          if (reset_moment) &
3396          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3397          count_reset_moment," steps"
3398          if (reset_vel) &
3399           write (iout,'(a,i10,a)') &
3400           "Velocities will be reset at random every",count_reset_vel,&
3401          " steps"
3402         endif
3403       else
3404         if(me.eq.king.or..not.out1file) &
3405          write (iout,'(a31)') "Microcanonical mode calculation"
3406       endif
3407       if(me.eq.king.or..not.out1file)then
3408        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3409        if (usampl) then
3410           write(iout,*) "MD running with constraints."
3411           write(iout,*) "Equilibration time ", eq_time, " mtus."
3412           write(iout,*) "Constraining ", nfrag," fragments."
3413           write(iout,*) "Length of each fragment, weight and q0:"
3414           do iset=1,nset
3415            write (iout,*) "Set of restraints #",iset
3416            do i=1,nfrag
3417               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3418                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3419            enddo
3420            write(iout,*) "constraints between ", npair, "fragments."
3421            write(iout,*) "constraint pairs, weights and q0:"
3422            do i=1,npair
3423             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3424                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3425            enddo
3426            write(iout,*) "angle constraints within ", nfrag_back,&
3427             "backbone fragments."
3428            write(iout,*) "fragment, weights:"
3429            do i=1,nfrag_back
3430             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3431                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3432                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3433            enddo
3434           enddo
3435         iset=mod(kolor,nset)+1
3436        endif
3437       endif
3438       if(me.eq.king.or..not.out1file) &
3439        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3440       return
3441       end subroutine read_MDpar
3442 !-----------------------------------------------------------------------------
3443       subroutine map_read
3444
3445       use map_data
3446 !      implicit real*8 (a-h,o-z)
3447 !      include 'DIMENSIONS'
3448 !      include 'COMMON.MAP'
3449 !      include 'COMMON.IOUNITS'
3450       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3451       character(len=80) :: mapcard      !,ucase
3452 !el local variables
3453       integer :: imap
3454 !     real(kind=8) :: var,ene
3455
3456       do imap=1,nmap
3457         read (inp,'(a)') mapcard
3458         mapcard=ucase(mapcard)
3459         if (index(mapcard,'PHI').gt.0) then
3460           kang(imap)=1
3461         else if (index(mapcard,'THE').gt.0) then
3462           kang(imap)=2
3463         else if (index(mapcard,'ALP').gt.0) then
3464           kang(imap)=3
3465         else if (index(mapcard,'OME').gt.0) then
3466           kang(imap)=4
3467         else
3468           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3469           stop 'Error - illegal variable spec in MAP card.'
3470         endif
3471         call readi (mapcard,'RES1',res1(imap),0)
3472         call readi (mapcard,'RES2',res2(imap),0)
3473         if (res1(imap).eq.0) then
3474           res1(imap)=res2(imap)
3475         else if (res2(imap).eq.0) then
3476           res2(imap)=res1(imap)
3477         endif
3478         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3479           write (iout,'(a)') &
3480           'Error - illegal definition of variable group in MAP.'
3481           stop 'Error - illegal definition of variable group in MAP.'
3482         endif
3483         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3484         call reada(mapcard,'TO',ang_to(imap),0.0D0)
3485         call readi(mapcard,'NSTEP',nstep(imap),0)
3486         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3487           write (iout,'(a)') &
3488            'Illegal boundary and/or step size specification in MAP.'
3489           stop 'Illegal boundary and/or step size specification in MAP.'
3490         endif
3491       enddo ! imap
3492       return
3493       end subroutine map_read
3494 !-----------------------------------------------------------------------------
3495       subroutine csaread
3496
3497       use control_data, only: vdisulf
3498       use csa_data
3499 !      implicit real*8 (a-h,o-z)
3500 !      include 'DIMENSIONS'
3501 !      include 'COMMON.IOUNITS'
3502 !      include 'COMMON.GEO'
3503 !      include 'COMMON.CSA'
3504 !      include 'COMMON.BANK'
3505 !      include 'COMMON.CONTROL'
3506 !      character(len=80) :: ucase
3507       character(len=620) :: mcmcard
3508 !el local variables
3509 !     integer :: ntf,ik,iw_pdb
3510 !     real(kind=8) :: var,ene
3511
3512       call card_concat(mcmcard,.true.)
3513
3514       call readi(mcmcard,'NCONF',nconf,50)
3515       call readi(mcmcard,'NADD',nadd,0)
3516       call readi(mcmcard,'JSTART',jstart,1)
3517       call readi(mcmcard,'JEND',jend,1)
3518       call readi(mcmcard,'NSTMAX',nstmax,500000)
3519       call readi(mcmcard,'N0',n0,1)
3520       call readi(mcmcard,'N1',n1,6)
3521       call readi(mcmcard,'N2',n2,4)
3522       call readi(mcmcard,'N3',n3,0)
3523       call readi(mcmcard,'N4',n4,0)
3524       call readi(mcmcard,'N5',n5,0)
3525       call readi(mcmcard,'N6',n6,10)
3526       call readi(mcmcard,'N7',n7,0)
3527       call readi(mcmcard,'N8',n8,0)
3528       call readi(mcmcard,'N9',n9,0)
3529       call readi(mcmcard,'N14',n14,0)
3530       call readi(mcmcard,'N15',n15,0)
3531       call readi(mcmcard,'N16',n16,0)
3532       call readi(mcmcard,'N17',n17,0)
3533       call readi(mcmcard,'N18',n18,0)
3534
3535       vdisulf=(index(mcmcard,'DYNSS').gt.0)
3536
3537       call readi(mcmcard,'NDIFF',ndiff,2)
3538       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3539       call readi(mcmcard,'IS1',is1,1)
3540       call readi(mcmcard,'IS2',is2,8)
3541       call readi(mcmcard,'NRAN0',nran0,4)
3542       call readi(mcmcard,'NRAN1',nran1,2)
3543       call readi(mcmcard,'IRR',irr,1)
3544       call readi(mcmcard,'NSEED',nseed,20)
3545       call readi(mcmcard,'NTOTAL',ntotal,10000)
3546       call reada(mcmcard,'CUT1',cut1,2.0d0)
3547       call reada(mcmcard,'CUT2',cut2,5.0d0)
3548       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3549       call readi(mcmcard,'ICMAX',icmax,3)
3550       call readi(mcmcard,'IRESTART',irestart,0)
3551 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
3552       ntbankm=0
3553 !!bankt
3554       call reada(mcmcard,'DELE',dele,20.0d0)
3555       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3556       call readi(mcmcard,'IREF',iref,0)
3557       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3558       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3559       call readi(mcmcard,'NCONF_IN',nconf_in,0)
3560       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3561       write (iout,*) "NCONF_IN",nconf_in
3562       return
3563       end subroutine csaread
3564 !-----------------------------------------------------------------------------
3565       subroutine mcmread
3566
3567       use mcm_data
3568       use control_data, only: MaxMoveType
3569       use MD_data
3570       use minim_data
3571 !      implicit real*8 (a-h,o-z)
3572 !      include 'DIMENSIONS'
3573 !      include 'COMMON.MCM'
3574 !      include 'COMMON.MCE'
3575 !      include 'COMMON.IOUNITS'
3576 !      character(len=80) :: ucase
3577       character(len=320) :: mcmcard
3578 !el local variables
3579       integer :: i
3580 !     real(kind=8) :: var,ene
3581
3582       call card_concat(mcmcard,.true.)
3583       call readi(mcmcard,'MAXACC',maxacc,100)
3584       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3585       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3586       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3587       call readi(mcmcard,'MAXREPM',maxrepm,200)
3588       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3589       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3590       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3591       call reada(mcmcard,'E_UP',e_up,5.0D0)
3592       call reada(mcmcard,'DELTE',delte,0.1D0)
3593       call readi(mcmcard,'NSWEEP',nsweep,5)
3594       call readi(mcmcard,'NSTEPH',nsteph,0)
3595       call readi(mcmcard,'NSTEPC',nstepc,0)
3596       call reada(mcmcard,'TMIN',tmin,298.0D0)
3597       call reada(mcmcard,'TMAX',tmax,298.0D0)
3598       call readi(mcmcard,'NWINDOW',nwindow,0)
3599       call readi(mcmcard,'PRINT_MC',print_mc,0)
3600       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3601       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3602       ent_read=(index(mcmcard,'ENT_READ').gt.0)
3603       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3604       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3605       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3606       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3607       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3608       if (nwindow.gt.0) then
3609         allocate(winstart(nwindow))     !!el (maxres)
3610         allocate(winend(nwindow))       !!el
3611         allocate(winlen(nwindow))       !!el
3612         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3613         do i=1,nwindow
3614           winlen(i)=winend(i)-winstart(i)+1
3615         enddo
3616       endif
3617       if (tmax.lt.tmin) tmax=tmin
3618       if (tmax.eq.tmin) then
3619         nstepc=0
3620         nsteph=0
3621       endif
3622       if (nstepc.gt.0 .and. nsteph.gt.0) then
3623         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3624         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3625       endif
3626       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3627 ! Probabilities of different move types
3628       sumpro_type(0)=0.0D0
3629       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3630       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3631       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3632       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
3633       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3634       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3635       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3636       do i=1,MaxMoveType
3637         print *,'i',i,' sumprotype',sumpro_type(i)
3638         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3639         print *,'i',i,' sumprotype',sumpro_type(i)
3640       enddo
3641       return
3642       end subroutine mcmread
3643 !-----------------------------------------------------------------------------
3644       subroutine read_minim
3645
3646       use minim_data
3647 !      implicit real*8 (a-h,o-z)
3648 !      include 'DIMENSIONS'
3649 !      include 'COMMON.MINIM'
3650 !      include 'COMMON.IOUNITS'
3651 !      character(len=80) :: ucase
3652       character(len=320) :: minimcard
3653 !el local variables
3654 !     integer :: ntf,ik,iw_pdb
3655 !     real(kind=8) :: var,ene
3656
3657       call card_concat(minimcard,.true.)
3658       call readi(minimcard,'MAXMIN',maxmin,2000)
3659       call readi(minimcard,'MAXFUN',maxfun,5000)
3660       call readi(minimcard,'MINMIN',minmin,maxmin)
3661       call readi(minimcard,'MINFUN',minfun,maxmin)
3662       call reada(minimcard,'TOLF',tolf,1.0D-2)
3663       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3664       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3665       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3666       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3667       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3668                'Options in energy minimization:'
3669       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3670        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3671        'MinMin:',MinMin,' MinFun:',MinFun,&
3672        ' TolF:',TolF,' RTolF:',RTolF
3673       return
3674       end subroutine read_minim
3675 !-----------------------------------------------------------------------------
3676       subroutine openunits
3677
3678       use MD_data, only: usampl
3679       use csa_data
3680       use MPI_data
3681       use control_data, only:out1file
3682       use control, only: getenv_loc
3683 !      implicit real*8 (a-h,o-z)
3684 !      include 'DIMENSIONS'
3685 #ifdef MPI
3686       include 'mpif.h'
3687       character(len=16) :: form,nodename
3688       integer :: nodelen,ierror,npos
3689 #endif
3690 !      include 'COMMON.SETUP'
3691 !      include 'COMMON.IOUNITS'
3692 !      include 'COMMON.MD'
3693 !      include 'COMMON.CONTROL'
3694       integer :: lenpre,lenpot,lentmp   !,ilen
3695 !el      external ilen
3696       character(len=3) :: out1file_text !,ucase
3697       character(len=3) :: ll
3698 !el      external ucase
3699 !el local variables
3700 !     integer :: ntf,ik,iw_pdb
3701 !     real(kind=8) :: var,ene
3702 !
3703 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3704       call getenv_loc("PREFIX",prefix)
3705       pref_orig = prefix
3706       call getenv_loc("POT",pot)
3707       call getenv_loc("DIRTMP",tmpdir)
3708       call getenv_loc("CURDIR",curdir)
3709       call getenv_loc("OUT1FILE",out1file_text)
3710 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3711       out1file_text=ucase(out1file_text)
3712       if (out1file_text(1:1).eq."Y") then
3713         out1file=.true.
3714       else
3715         out1file=fg_rank.gt.0
3716       endif
3717       lenpre=ilen(prefix)
3718       lenpot=ilen(pot)
3719       lentmp=ilen(tmpdir)
3720       if (lentmp.gt.0) then
3721           write (*,'(80(1h!))')
3722           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
3723           write (*,'(80(1h!))')
3724           write (*,*)"All output files will be on node /tmp directory."
3725 #ifdef MPI
3726         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3727         if (me.eq.king) then
3728           write (*,*) "The master node is ",nodename
3729         else if (fg_rank.eq.0) then
3730           write (*,*) "I am the CG slave node ",nodename
3731         else
3732           write (*,*) "I am the FG slave node ",nodename
3733         endif
3734 #endif
3735         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3736         lenpre = lentmp+lenpre+1
3737       endif
3738       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3739 ! Get the names and open the input files
3740 #if defined(WINIFL) || defined(WINPGI)
3741       open(1,file=pref_orig(:ilen(pref_orig))// &
3742         '.inp',status='old',readonly,shared)
3743        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3744 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3745 ! Get parameter filenames and open the parameter files.
3746       call getenv_loc('BONDPAR',bondname)
3747       open (ibond,file=bondname,status='old',readonly,shared)
3748       call getenv_loc('THETPAR',thetname)
3749       open (ithep,file=thetname,status='old',readonly,shared)
3750       call getenv_loc('ROTPAR',rotname)
3751       open (irotam,file=rotname,status='old',readonly,shared)
3752       call getenv_loc('TORPAR',torname)
3753       open (itorp,file=torname,status='old',readonly,shared)
3754       call getenv_loc('TORDPAR',tordname)
3755       open (itordp,file=tordname,status='old',readonly,shared)
3756       call getenv_loc('FOURIER',fouriername)
3757       open (ifourier,file=fouriername,status='old',readonly,shared)
3758       call getenv_loc('ELEPAR',elename)
3759       open (ielep,file=elename,status='old',readonly,shared)
3760       call getenv_loc('SIDEPAR',sidename)
3761       open (isidep,file=sidename,status='old',readonly,shared)
3762 #elif (defined CRAY) || (defined AIX)
3763       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3764         action='read')
3765 !      print *,"Processor",myrank," opened file 1"
3766       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3767 !      print *,"Processor",myrank," opened file 9"
3768 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3769 ! Get parameter filenames and open the parameter files.
3770       call getenv_loc('BONDPAR',bondname)
3771       open (ibond,file=bondname,status='old',action='read')
3772 !      print *,"Processor",myrank," opened file IBOND"
3773       call getenv_loc('THETPAR',thetname)
3774       open (ithep,file=thetname,status='old',action='read')
3775 !      print *,"Processor",myrank," opened file ITHEP"
3776       call getenv_loc('ROTPAR',rotname)
3777       open (irotam,file=rotname,status='old',action='read')
3778 !      print *,"Processor",myrank," opened file IROTAM"
3779       call getenv_loc('TORPAR',torname)
3780       open (itorp,file=torname,status='old',action='read')
3781 !      print *,"Processor",myrank," opened file ITORP"
3782       call getenv_loc('TORDPAR',tordname)
3783       open (itordp,file=tordname,status='old',action='read')
3784 !      print *,"Processor",myrank," opened file ITORDP"
3785       call getenv_loc('SCCORPAR',sccorname)
3786       open (isccor,file=sccorname,status='old',action='read')
3787 !      print *,"Processor",myrank," opened file ISCCOR"
3788       call getenv_loc('FOURIER',fouriername)
3789       open (ifourier,file=fouriername,status='old',action='read')
3790 !      print *,"Processor",myrank," opened file IFOURIER"
3791       call getenv_loc('ELEPAR',elename)
3792       open (ielep,file=elename,status='old',action='read')
3793 !      print *,"Processor",myrank," opened file IELEP"
3794       call getenv_loc('SIDEPAR',sidename)
3795       open (isidep,file=sidename,status='old',action='read')
3796       call getenv_loc('LIPTRANPAR',liptranname)
3797       open (iliptranpar,file=liptranname,status='old',action='read')
3798       call getenv_loc('TUBEPAR',tubename)
3799       open (itube,file=tubename,status='old',action='read')
3800
3801 !      print *,"Processor",myrank," opened file ISIDEP"
3802 !      print *,"Processor",myrank," opened parameter files"
3803 #elif (defined G77)
3804       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3805       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3806 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3807 ! Get parameter filenames and open the parameter files.
3808       call getenv_loc('BONDPAR',bondname)
3809       open (ibond,file=bondname,status='old')
3810       call getenv_loc('THETPAR',thetname)
3811       open (ithep,file=thetname,status='old')
3812       call getenv_loc('ROTPAR',rotname)
3813       open (irotam,file=rotname,status='old')
3814       call getenv_loc('TORPAR',torname)
3815       open (itorp,file=torname,status='old')
3816       call getenv_loc('TORDPAR',tordname)
3817       open (itordp,file=tordname,status='old')
3818       call getenv_loc('SCCORPAR',sccorname)
3819       open (isccor,file=sccorname,status='old')
3820       call getenv_loc('FOURIER',fouriername)
3821       open (ifourier,file=fouriername,status='old')
3822       call getenv_loc('ELEPAR',elename)
3823       open (ielep,file=elename,status='old')
3824       call getenv_loc('SIDEPAR',sidename)
3825       open (isidep,file=sidename,status='old')
3826       call getenv_loc('LIPTRANPAR',liptranname)
3827       open (iliptranpar,file=liptranname,status='old')
3828       call getenv_loc('TUBEPAR',tubename)
3829       open (itube,file=tubename,status='old')
3830 #else
3831       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3832         readonly)
3833        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3834 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3835 ! Get parameter filenames and open the parameter files.
3836       call getenv_loc('BONDPAR',bondname)
3837       open (ibond,file=bondname,status='old',action='read')
3838       call getenv_loc('THETPAR',thetname)
3839       open (ithep,file=thetname,status='old',action='read')
3840       call getenv_loc('ROTPAR',rotname)
3841       open (irotam,file=rotname,status='old',action='read')
3842       call getenv_loc('TORPAR',torname)
3843       open (itorp,file=torname,status='old',action='read')
3844       call getenv_loc('TORDPAR',tordname)
3845       open (itordp,file=tordname,status='old',action='read')
3846       call getenv_loc('SCCORPAR',sccorname)
3847       open (isccor,file=sccorname,status='old',action='read')
3848 #ifndef CRYST_THETA
3849       call getenv_loc('THETPARPDB',thetname_pdb)
3850       print *,"thetname_pdb ",thetname_pdb
3851       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3852       print *,ithep_pdb," opened"
3853 #endif
3854       call getenv_loc('FOURIER',fouriername)
3855       open (ifourier,file=fouriername,status='old',readonly)
3856       call getenv_loc('ELEPAR',elename)
3857       open (ielep,file=elename,status='old',readonly)
3858       call getenv_loc('SIDEPAR',sidename)
3859       open (isidep,file=sidename,status='old',readonly)
3860       call getenv_loc('LIPTRANPAR',liptranname)
3861       open (iliptranpar,file=liptranname,status='old',action='read')
3862       call getenv_loc('TUBEPAR',tubename)
3863       open (itube,file=tubename,status='old',action='read')
3864
3865 #ifndef CRYST_SC
3866       call getenv_loc('ROTPARPDB',rotname_pdb)
3867       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3868 #endif
3869 #endif
3870 #ifndef OLDSCP
3871 !
3872 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3873 ! Use -DOLDSCP to use hard-coded constants instead.
3874 !
3875       call getenv_loc('SCPPAR',scpname)
3876 #if defined(WINIFL) || defined(WINPGI)
3877       open (iscpp,file=scpname,status='old',readonly,shared)
3878 #elif (defined CRAY)  || (defined AIX)
3879       open (iscpp,file=scpname,status='old',action='read')
3880 #elif (defined G77)
3881       open (iscpp,file=scpname,status='old')
3882 #else
3883       open (iscpp,file=scpname,status='old',action='read')
3884 #endif
3885 #endif
3886       call getenv_loc('PATTERN',patname)
3887 #if defined(WINIFL) || defined(WINPGI)
3888       open (icbase,file=patname,status='old',readonly,shared)
3889 #elif (defined CRAY)  || (defined AIX)
3890       open (icbase,file=patname,status='old',action='read')
3891 #elif (defined G77)
3892       open (icbase,file=patname,status='old')
3893 #else
3894       open (icbase,file=patname,status='old',action='read')
3895 #endif
3896 #ifdef MPI
3897 ! Open output file only for CG processes
3898 !      print *,"Processor",myrank," fg_rank",fg_rank
3899       if (fg_rank.eq.0) then
3900
3901       if (nodes.eq.1) then
3902         npos=3
3903       else
3904         npos = dlog10(dfloat(nodes-1))+1
3905       endif
3906       if (npos.lt.3) npos=3
3907       write (liczba,'(i1)') npos
3908       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3909         //')'
3910       write (liczba,form) me
3911       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3912         liczba(:ilen(liczba))
3913       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3914         //'.int'
3915       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3916         //'.pdb'
3917       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3918         liczba(:ilen(liczba))//'.mol2'
3919       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3920         liczba(:ilen(liczba))//'.stat'
3921       if (lentmp.gt.0) &
3922         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3923             //liczba(:ilen(liczba))//'.stat')
3924       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3925         //'.rst'
3926       if(usampl) then
3927           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3928        liczba(:ilen(liczba))//'.const'
3929       endif
3930
3931       endif
3932 #else
3933       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3934       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3935       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3936       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3937       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3938       if (lentmp.gt.0) &
3939         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3940           '.stat')
3941       rest2name=prefix(:ilen(prefix))//'.rst'
3942       if(usampl) then
3943          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3944       endif
3945 #endif
3946 #if defined(AIX) || defined(PGI)
3947       if (me.eq.king .or. .not. out1file) &
3948          open(iout,file=outname,status='unknown')
3949 #ifdef DEBUG
3950       if (fg_rank.gt.0) then
3951         write (liczba,'(i3.3)') myrank/nfgtasks
3952         write (ll,'(bz,i3.3)') fg_rank
3953         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3954          status='unknown')
3955       endif
3956 #endif
3957       if(me.eq.king) then
3958        open(igeom,file=intname,status='unknown',position='append')
3959        open(ipdb,file=pdbname,status='unknown')
3960        open(imol2,file=mol2name,status='unknown')
3961        open(istat,file=statname,status='unknown',position='append')
3962       else
3963 !1out       open(iout,file=outname,status='unknown')
3964       endif
3965 #else
3966       if (me.eq.king .or. .not.out1file) &
3967           open(iout,file=outname,status='unknown')
3968 #ifdef DEBUG
3969       if (fg_rank.gt.0) then
3970         write (liczba,'(i3.3)') myrank/nfgtasks
3971         write (ll,'(bz,i3.3)') fg_rank
3972         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3973          status='unknown')
3974       endif
3975 #endif
3976       if(me.eq.king) then
3977        open(igeom,file=intname,status='unknown',access='append')
3978        open(ipdb,file=pdbname,status='unknown')
3979        open(imol2,file=mol2name,status='unknown')
3980        open(istat,file=statname,status='unknown',access='append')
3981       else
3982 !1out       open(iout,file=outname,status='unknown')
3983       endif
3984 #endif
3985       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3986       csa_seed=prefix(:lenpre)//'.CSA.seed'
3987       csa_history=prefix(:lenpre)//'.CSA.history'
3988       csa_bank=prefix(:lenpre)//'.CSA.bank'
3989       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3990       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3991       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3992 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3993       csa_int=prefix(:lenpre)//'.int'
3994       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3995       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3996       csa_in=prefix(:lenpre)//'.CSA.in'
3997 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3998 ! Write file names
3999       if (me.eq.king)then
4000       write (iout,'(80(1h-))')
4001       write (iout,'(30x,a)') "FILE ASSIGNMENT"
4002       write (iout,'(80(1h-))')
4003       write (iout,*) "Input file                      : ",&
4004         pref_orig(:ilen(pref_orig))//'.inp'
4005       write (iout,*) "Output file                     : ",&
4006         outname(:ilen(outname))
4007       write (iout,*)
4008       write (iout,*) "Sidechain potential file        : ",&
4009         sidename(:ilen(sidename))
4010 #ifndef OLDSCP
4011       write (iout,*) "SCp potential file              : ",&
4012         scpname(:ilen(scpname))
4013 #endif
4014       write (iout,*) "Electrostatic potential file    : ",&
4015         elename(:ilen(elename))
4016       write (iout,*) "Cumulant coefficient file       : ",&
4017         fouriername(:ilen(fouriername))
4018       write (iout,*) "Torsional parameter file        : ",&
4019         torname(:ilen(torname))
4020       write (iout,*) "Double torsional parameter file : ",&
4021         tordname(:ilen(tordname))
4022       write (iout,*) "SCCOR parameter file : ",&
4023         sccorname(:ilen(sccorname))
4024       write (iout,*) "Bond & inertia constant file    : ",&
4025         bondname(:ilen(bondname))
4026       write (iout,*) "Bending parameter file          : ",&
4027         thetname(:ilen(thetname))
4028       write (iout,*) "Rotamer parameter file          : ",&
4029         rotname(:ilen(rotname))
4030 !el----
4031 #ifndef CRYST_THETA
4032       write (iout,*) "Thetpdb parameter file          : ",&
4033         thetname_pdb(:ilen(thetname_pdb))
4034 #endif
4035 !el
4036       write (iout,*) "Threading database              : ",&
4037         patname(:ilen(patname))
4038       if (lentmp.ne.0) &
4039       write (iout,*)" DIRTMP                          : ",&
4040         tmpdir(:lentmp)
4041       write (iout,'(80(1h-))')
4042       endif
4043       return
4044       end subroutine openunits
4045 !-----------------------------------------------------------------------------
4046       subroutine readrst
4047
4048       use geometry_data, only: nres,dc
4049       use MD_data
4050 !      implicit real*8 (a-h,o-z)
4051 !      include 'DIMENSIONS'
4052 !      include 'COMMON.CHAIN'
4053 !      include 'COMMON.IOUNITS'
4054 !      include 'COMMON.MD'
4055 !el local variables
4056       integer ::i,j
4057 !     real(kind=8) :: var,ene
4058
4059       open(irest2,file=rest2name,status='unknown')
4060       read(irest2,*) totT,EK,potE,totE,t_bath
4061       totTafm=totT
4062 !      do i=1,2*nres
4063 ! AL 4/17/17: Now reading d_t(0,:) too
4064       do i=0,2*nres
4065          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4066       enddo
4067 !      do i=1,2*nres
4068 ! AL 4/17/17: Now reading d_c(0,:) too
4069       do i=0,2*nres
4070          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4071       enddo
4072       if(usampl) then
4073              read (irest2,*) iset
4074       endif
4075       close(irest2)
4076       return
4077       end subroutine readrst
4078 !-----------------------------------------------------------------------------
4079       subroutine read_fragments
4080
4081       use energy_data
4082 !      use geometry
4083       use control_data, only:out1file
4084       use MD_data
4085       use MPI_data
4086 !      implicit real*8 (a-h,o-z)
4087 !      include 'DIMENSIONS'
4088 #ifdef MPI
4089       include 'mpif.h'
4090 #endif
4091 !      include 'COMMON.SETUP'
4092 !      include 'COMMON.CHAIN'
4093 !      include 'COMMON.IOUNITS'
4094 !      include 'COMMON.MD'
4095 !      include 'COMMON.CONTROL'
4096 !el local variables
4097       integer :: i
4098 !     real(kind=8) :: var,ene
4099
4100       read(inp,*) nset,nfrag,npair,nfrag_back
4101
4102 !el from module energy
4103 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
4104       if(.not.allocated(wfrag_back)) then
4105         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4106         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4107
4108         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4109         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
4110
4111         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4112         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
4113       endif
4114
4115       if(me.eq.king.or..not.out1file) &
4116        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4117         " nfrag_back",nfrag_back
4118       do iset=1,nset
4119          read(inp,*) mset(iset)
4120        do i=1,nfrag
4121          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4122            qinfrag(i,iset)
4123          if(me.eq.king.or..not.out1file) &
4124           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4125            ifrag(2,i,iset), qinfrag(i,iset)
4126        enddo
4127        do i=1,npair
4128         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4129           qinpair(i,iset)
4130         if(me.eq.king.or..not.out1file) &
4131          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4132           ipair(2,i,iset), qinpair(i,iset)
4133        enddo
4134        do i=1,nfrag_back
4135         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4136            wfrag_back(3,i,iset),&
4137            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4138         if(me.eq.king.or..not.out1file) &
4139          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4140          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4141        enddo
4142       enddo
4143       return
4144       end subroutine read_fragments
4145 !-----------------------------------------------------------------------------
4146 ! shift.F       io_csa
4147 !-----------------------------------------------------------------------------
4148       subroutine csa_read
4149
4150       use csa_data
4151 !      implicit real*8 (a-h,o-z)
4152 !      include 'DIMENSIONS'
4153 !      include 'COMMON.CSA'
4154 !      include 'COMMON.BANK'
4155 !      include 'COMMON.IOUNITS'
4156 !el local variables
4157 !     integer :: ntf,ik,iw_pdb
4158 !     real(kind=8) :: var,ene
4159
4160       open(icsa_in,file=csa_in,status="old",err=100)
4161        read(icsa_in,*) nconf
4162        read(icsa_in,*) jstart,jend
4163        read(icsa_in,*) nstmax
4164        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4165        read(icsa_in,*) nran0,nran1,irr
4166        read(icsa_in,*) nseed
4167        read(icsa_in,*) ntotal,cut1,cut2
4168        read(icsa_in,*) estop
4169        read(icsa_in,*) icmax,irestart
4170        read(icsa_in,*) ntbankm,dele,difcut
4171        read(icsa_in,*) iref,rmscut,pnccut
4172        read(icsa_in,*) ndiff
4173       close(icsa_in)
4174
4175       return
4176
4177  100  continue
4178       return
4179       end subroutine csa_read
4180 !-----------------------------------------------------------------------------
4181       subroutine initial_write
4182
4183       use csa_data
4184 !      implicit real*8 (a-h,o-z)
4185 !      include 'DIMENSIONS'
4186 !      include 'COMMON.CSA'
4187 !      include 'COMMON.BANK'
4188 !      include 'COMMON.IOUNITS'
4189 !el local variables
4190 !     integer :: ntf,ik,iw_pdb
4191 !     real(kind=8) :: var,ene
4192
4193       open(icsa_seed,file=csa_seed,status="unknown")
4194        write(icsa_seed,*) "seed"
4195       close(31)
4196 #if defined(AIX) || defined(PGI)
4197        open(icsa_history,file=csa_history,status="unknown",&
4198         position="append")
4199 #else
4200        open(icsa_history,file=csa_history,status="unknown",&
4201         access="append")
4202 #endif
4203        write(icsa_history,*) nconf
4204        write(icsa_history,*) jstart,jend
4205        write(icsa_history,*) nstmax
4206        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4207        write(icsa_history,*) nran0,nran1,irr
4208        write(icsa_history,*) nseed
4209        write(icsa_history,*) ntotal,cut1,cut2
4210        write(icsa_history,*) estop
4211        write(icsa_history,*) icmax,irestart
4212        write(icsa_history,*) ntbankm,dele,difcut
4213        write(icsa_history,*) iref,rmscut,pnccut
4214        write(icsa_history,*) ndiff
4215
4216        write(icsa_history,*)
4217       close(icsa_history)
4218
4219       open(icsa_bank1,file=csa_bank1,status="unknown")
4220        write(icsa_bank1,*) 0
4221       close(icsa_bank1)
4222
4223       return
4224       end subroutine initial_write
4225 !-----------------------------------------------------------------------------
4226       subroutine restart_write
4227
4228       use csa_data
4229 !      implicit real*8 (a-h,o-z)
4230 !      include 'DIMENSIONS'
4231 !      include 'COMMON.IOUNITS'
4232 !      include 'COMMON.CSA'
4233 !      include 'COMMON.BANK'
4234 !el local variables
4235 !     integer :: ntf,ik,iw_pdb
4236 !     real(kind=8) :: var,ene
4237
4238 #if defined(AIX) || defined(PGI)
4239        open(icsa_history,file=csa_history,position="append")
4240 #else
4241        open(icsa_history,file=csa_history,access="append")
4242 #endif
4243        write(icsa_history,*)
4244        write(icsa_history,*) "This is restart"
4245        write(icsa_history,*)
4246        write(icsa_history,*) nconf
4247        write(icsa_history,*) jstart,jend
4248        write(icsa_history,*) nstmax
4249        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4250        write(icsa_history,*) nran0,nran1,irr
4251        write(icsa_history,*) nseed
4252        write(icsa_history,*) ntotal,cut1,cut2
4253        write(icsa_history,*) estop
4254        write(icsa_history,*) icmax,irestart
4255        write(icsa_history,*) ntbankm,dele,difcut
4256        write(icsa_history,*) iref,rmscut,pnccut
4257        write(icsa_history,*) ndiff
4258        write(icsa_history,*)
4259        write(icsa_history,*) "irestart is: ", irestart
4260
4261        write(icsa_history,*)
4262       close(icsa_history)
4263
4264       return
4265       end subroutine restart_write
4266 !-----------------------------------------------------------------------------
4267 ! test.F
4268 !-----------------------------------------------------------------------------
4269       subroutine write_pdb(npdb,titelloc,ee)
4270
4271 !      implicit real*8 (a-h,o-z)
4272 !      include 'DIMENSIONS'
4273 !      include 'COMMON.IOUNITS'
4274       character(len=50) :: titelloc1
4275       character*(*) :: titelloc
4276       character(len=3) :: zahl
4277       character(len=5) :: liczba5
4278       real(kind=8) :: ee
4279       integer :: npdb   !,ilen
4280 !el      external ilen
4281 !el local variables
4282       integer :: lenpre
4283 !     real(kind=8) :: var,ene
4284
4285       titelloc1=titelloc
4286       lenpre=ilen(prefix)
4287       if (npdb.lt.1000) then
4288        call numstr(npdb,zahl)
4289        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4290       else
4291         if (npdb.lt.10000) then
4292          write(liczba5,'(i1,i4)') 0,npdb
4293         else
4294          write(liczba5,'(i5)') npdb
4295         endif
4296         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4297       endif
4298       call pdbout(ee,titelloc1,ipdb)
4299       close(ipdb)
4300       return
4301       end subroutine write_pdb
4302 !-----------------------------------------------------------------------------
4303 ! thread.F
4304 !-----------------------------------------------------------------------------
4305       subroutine write_thread_summary
4306 ! Thread the sequence through a database of known structures
4307       use control_data, only: refstr
4308 !      use geometry
4309       use energy_data, only: n_ene_comp
4310       use compare_data
4311 !      implicit real*8 (a-h,o-z)
4312 !      include 'DIMENSIONS'
4313 #ifdef MPI
4314       use MPI_data      !include 'COMMON.INFO'
4315       include 'mpif.h'
4316 #endif
4317 !      include 'COMMON.CONTROL'
4318 !      include 'COMMON.CHAIN'
4319 !      include 'COMMON.DBASE'
4320 !      include 'COMMON.INTERACT'
4321 !      include 'COMMON.VAR'
4322 !      include 'COMMON.THREAD'
4323 !      include 'COMMON.FFIELD'
4324 !      include 'COMMON.SBRIDGE'
4325 !      include 'COMMON.HEADER'
4326 !      include 'COMMON.NAMES'
4327 !      include 'COMMON.IOUNITS'
4328 !      include 'COMMON.TIME1'
4329
4330       integer,dimension(maxthread) :: ip
4331       real(kind=8),dimension(0:n_ene) :: energia
4332 !el local variables
4333       integer :: i,j,ii,jj,ipj,ik,kk,ist
4334       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4335
4336       write (iout,'(30x,a/)') &
4337        '  *********** Summary threading statistics ************'
4338       write (iout,'(a)') 'Initial energies:'
4339       write (iout,'(a4,2x,a12,14a14,3a8)') &
4340        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4341        'RMSnat','NatCONT','NNCONT','RMS'
4342 ! Energy sort patterns
4343       do i=1,nthread
4344         ip(i)=i
4345       enddo
4346       do i=1,nthread-1
4347         enet=ener(n_ene-1,ip(i))
4348         jj=i
4349         do j=i+1,nthread
4350           if (ener(n_ene-1,ip(j)).lt.enet) then
4351             jj=j
4352             enet=ener(n_ene-1,ip(j))
4353           endif
4354         enddo
4355         if (jj.ne.i) then
4356           ipj=ip(jj)
4357           ip(jj)=ip(i)
4358           ip(i)=ipj
4359         endif
4360       enddo
4361       do ik=1,nthread
4362         i=ip(ik)
4363         ii=ipatt(1,i)
4364         ist=nres_base(2,ii)+ipatt(2,i)
4365         do kk=1,n_ene_comp
4366           energia(i)=ener0(kk,i)
4367         enddo
4368         etot=ener0(n_ene_comp+1,i)
4369         rmsnat=ener0(n_ene_comp+2,i)
4370         rms=ener0(n_ene_comp+3,i)
4371         frac=ener0(n_ene_comp+4,i)
4372         frac_nn=ener0(n_ene_comp+5,i)
4373
4374         if (refstr) then
4375         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4376         i,str_nam(ii),ist+1,&
4377         (energia(print_order(kk)),kk=1,nprint_ene),&
4378         etot,rmsnat,frac,frac_nn,rms
4379         else
4380         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4381         i,str_nam(ii),ist+1,&
4382         (energia(print_order(kk)),kk=1,nprint_ene),etot
4383         endif
4384       enddo
4385       write (iout,'(//a)') 'Final energies:'
4386       write (iout,'(a4,2x,a12,17a14,3a8)') &
4387        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4388        'RMSnat','NatCONT','NNCONT','RMS'
4389       do ik=1,nthread
4390         i=ip(ik)
4391         ii=ipatt(1,i)
4392         ist=nres_base(2,ii)+ipatt(2,i)
4393         do kk=1,n_ene_comp
4394           energia(kk)=ener(kk,ik)
4395         enddo
4396         etot=ener(n_ene_comp+1,i)
4397         rmsnat=ener(n_ene_comp+2,i)
4398         rms=ener(n_ene_comp+3,i)
4399         frac=ener(n_ene_comp+4,i)
4400         frac_nn=ener(n_ene_comp+5,i)
4401         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4402         i,str_nam(ii),ist+1,&
4403         (energia(print_order(kk)),kk=1,nprint_ene),&
4404         etot,rmsnat,frac,frac_nn,rms
4405       enddo
4406       write (iout,'(/a/)') 'IEXAM array:'
4407       write (iout,'(i5)') nexcl
4408       do i=1,nexcl
4409         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4410       enddo
4411       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4412        'Max. time for threading step ',max_time_for_thread,&
4413        'Average time for threading step: ',ave_time_for_thread
4414       return
4415       end subroutine write_thread_summary
4416 !-----------------------------------------------------------------------------
4417       subroutine write_stat_thread(ithread,ipattern,ist)
4418
4419       use energy_data, only: n_ene_comp
4420       use compare_data
4421 !      implicit real*8 (a-h,o-z)
4422 !      include "DIMENSIONS"
4423 !      include "COMMON.CONTROL"
4424 !      include "COMMON.IOUNITS"
4425 !      include "COMMON.THREAD"
4426 !      include "COMMON.FFIELD"
4427 !      include "COMMON.DBASE"
4428 !      include "COMMON.NAMES"
4429       real(kind=8),dimension(0:n_ene) :: energia
4430 !el local variables
4431       integer :: ithread,ipattern,ist,i
4432       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4433
4434 #if defined(AIX) || defined(PGI)
4435       open(istat,file=statname,position='append')
4436 #else
4437       open(istat,file=statname,access='append')
4438 #endif
4439       do i=1,n_ene_comp
4440         energia(i)=ener(i,ithread)
4441       enddo
4442       etot=ener(n_ene_comp+1,ithread)
4443       rmsnat=ener(n_ene_comp+2,ithread)
4444       rms=ener(n_ene_comp+3,ithread)
4445       frac=ener(n_ene_comp+4,ithread)
4446       frac_nn=ener(n_ene_comp+5,ithread)
4447       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4448         ithread,str_nam(ipattern),ist+1,&
4449         (energia(print_order(i)),i=1,nprint_ene),&
4450         etot,rmsnat,frac,frac_nn,rms
4451       close (istat)
4452       return
4453       end subroutine write_stat_thread
4454 !-----------------------------------------------------------------------------
4455 #endif
4456 !-----------------------------------------------------------------------------
4457       end module io_config