source/unres/io_config.f90

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