source/unres/io_config.f90

   1       module io_config
   2
   3       use names
   4       use io_units
   5       use io_base
   6       use geometry_data
   7       use geometry
   8       use control_data, only:maxterm_sccor
   9       implicit none
  10 !-----------------------------------------------------------------------------
  11 ! Max. number of residue types and parameters in expressions for
  12 ! virtual-bond angle bending potentials
  13 !      integer,parameter :: maxthetyp=3
  14 !      integer,parameter :: maxthetyp1=maxthetyp+1
  15 !     ,maxtheterm=20,
  16 !     & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
  17 !     & mmaxtheterm=maxtheterm)
  18 !-----------------------------------------------------------------------------
  19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
  20 ! and the number of terms in double torsionals
  21 !      integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
  22 !      parameter (maxtor=4,maxterm=10)
  23 !-----------------------------------------------------------------------------
  24 ! Max number of torsional terms in SCCOR
  25 !el      integer,parameter :: maxterm_sccor=6
  26 !-----------------------------------------------------------------------------
  27       character(len=1),dimension(:),allocatable :: secstruc     !(maxres)
  28 !-----------------------------------------------------------------------------
  29 !
  30 !
  31 !-----------------------------------------------------------------------------
  32       contains
  33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
  34 !-----------------------------------------------------------------------------
  35 ! bank.F    io_csa
  36 !-----------------------------------------------------------------------------
  37       subroutine write_rbank(jlee,adif,nft)
  38
  39       use csa_data
  40       use geometry_data, only: nres,rad2deg
  41 !      implicit real*8 (a-h,o-z)
  42 !      include 'DIMENSIONS'
  43 !      include 'COMMON.IOUNITS'
  44 !      include 'COMMON.CSA'
  45 !      include 'COMMON.BANK'
  46 !      include 'COMMON.CHAIN'
  47 !      include 'COMMON.GEO'
  48 !el local variables
  49       integer :: nft,i,k,j,l,jlee
  50       real(kind=8) :: adif
  51
  52       open(icsa_rbank,file=csa_rbank,status="unknown")
  53       write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
  54       do k=1,nbank
  55        write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
  56        do j=1,numch
  57         do l=2,nres-1
  58          write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
  59         enddo
  60        enddo
  61       enddo
  62       close(icsa_rbank)
  63
  64   850 format (10f8.3)
  65   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
  66               i8,i10,i2,f15.5)
  67   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
  68                   ' %NC ',0pf5.2)
  69
  70       return
  71       end subroutine write_rbank
  72 !-----------------------------------------------------------------------------
  73       subroutine read_rbank(jlee,adif)
  74
  75       use csa_data
  76       use geometry_data, only: nres,deg2rad
  77       use MPI_data
  78 !      implicit real*8 (a-h,o-z)
  79 !      include 'DIMENSIONS'
  80       include 'mpif.h'
  81 !      include 'COMMON.IOUNITS'
  82 !      include 'COMMON.CSA'
  83 !      include 'COMMON.BANK'
  84 !      include 'COMMON.CHAIN'
  85 !      include 'COMMON.GEO'
  86 !      include 'COMMON.SETUP'
  87       character(len=80) :: karta
  88 !el local variables
  89       integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
  90                  ierror,ierrcode,jlee,jleer
  91       real(kind=8) :: adif
  92
  93       open(icsa_rbank,file=csa_rbank,status="old")
  94       read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
  95       print *,jleer,nbankr,nstepr,nftr,icycler,adif
  96 !       print *, 'adif from read_rbank ',adif
  97 #ifdef MPI
  98       if(nbankr.ne.nbank) then
  99         write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
 100         ' NBANK',nbank
 101         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 102       endif
 103       if(jleer.ne.jlee) then
 104         write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
 105         ' JLEE',jlee
 106         call mpi_abort(mpi_comm_world,ierror,ierrcode)
 107       endif
 108 #endif
 109
 110       kk=0
 111       do k=1,nbankr
 112         read (icsa_rbank,'(a80)') karta
 113         write(iout,*) "READ_RBANK: kk=",kk
 114         write(iout,*) karta
 115 !        if (index(karta,"*").gt.0) then
 116 !          write (iout,*) "***** Stars in bankr ***** k=",k,
 117 !     &      " skipped"
 118 !          do j=1,numch
 119 !            do l=2,nres-1
 120 !              read (30,850) (rdummy,i=1,4)
 121 !            enddo
 122 !          enddo
 123 !        else
 124           kk=kk+1
 125           call reada(karta,"total E",rene(kk),1.0d20)
 126           call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
 127           call reada(karta,"%NC",rpncn(kk),0.0d0)
 128           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
 129             "%NC",bpncn(kk),ibank(kk)
 130 !          read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
 131           do j=1,numch
 132             do l=2,nres-1
 133               read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
 134 !              write (iout,850) (rvar(i,l,j,kk),i=1,4)
 135               do i=1,4
 136                 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
 137               enddo
 138             enddo
 139           enddo
 140 !        endif
 141       enddo
 142 !d      write (*,*) "read_rbank ******************* kk",kk,
 143 !d     &  "nbankr",nbankr
 144       if (kk.lt.nbankr) nbankr=kk
 145 !d      do kk=1,nbankr
 146 !d          print *,"kk=",kk
 147 !d          do j=1,numch
 148 !d            do l=2,nres-1
 149 !d              write (*,850) (rvar(i,l,j,kk),i=1,4)
 150 !d            enddo
 151 !d          enddo
 152 !d      enddo
 153       close(icsa_rbank)
 154
 155   850 format (10f8.3)
 156   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 157   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
 158
 159       return
 160       end subroutine read_rbank
 161 !-----------------------------------------------------------------------------
 162       subroutine write_bank(jlee,nft)
 163
 164       use csa_data
 165       use control_data, only: vdisulf
 166       use geometry_data, only: nres,rad2deg
 167 !      implicit real*8 (a-h,o-z)
 168 !      include 'DIMENSIONS'
 169 !      include 'COMMON.IOUNITS'
 170 !      include 'COMMON.CSA'
 171 !      include 'COMMON.BANK'
 172 !      include 'COMMON.CHAIN'
 173 !      include 'COMMON.GEO'
 174 !      include 'COMMON.SBRIDGE'
 175 !     include 'COMMON.CONTROL'
 176       character(len=7) :: chtmp
 177       character(len=40) :: chfrm
 178 !el      external ilen
 179 !el local variables
 180       integer :: nft,k,l,i,j,jlee
 181
 182       open(icsa_bank,file=csa_bank,status="unknown")
 183       write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
 184       write (icsa_bank,902) nglob_csa, eglob_csa
 185       open (igeom,file=intname,status='UNKNOWN')
 186       do k=1,nbank
 187        write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 188        if (vdisulf) write (icsa_bank,'(101i4)') &
 189           bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
 190        do j=1,numch
 191         do l=2,nres-1
 192          write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 193         enddo
 194        enddo
 195        if (bvar_nss(k).le.9) then
 196          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 197            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 198        else
 199          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 200            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 201          write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
 202                                       bvar_ss(2,i,k),i=10,bvar_nss(k))
 203        endif
 204        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 205        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 206        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 207        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 208       enddo
 209       close(icsa_bank)
 210       close(igeom)
 211
 212       if (nstep/200.gt.ilastnstep) then
 213
 214        ilastnstep=(ilastnstep+1)*1.5
 215        write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
 216        write(chtmp,chfrm) nstep
 217        open(icsa_int,file=prefix(:ilen(prefix)) &
 218                //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
 219        do k=1,nbank
 220         if (bvar_nss(k).le.9) then
 221          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 222         bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
 223         else
 224          write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 225            bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
 226          write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
 227                                 bvar_ss(2,i,k),i=10,bvar_nss(k))
 228         endif
 229         write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 230         write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 231         write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 232         write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 233        enddo
 234        close(icsa_int)
 235       endif
 236
 237
 238   200 format (8f10.4)
 239   850 format (10f8.3)
 240   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
 241               i8,i10,i2,f15.5)
 242   902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
 243   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 244               ' %NC ',0pf5.2,i5)
 245
 246       return
 247       end subroutine write_bank
 248 !-----------------------------------------------------------------------------
 249       subroutine write_bank_reminimized(jlee,nft)
 250
 251       use csa_data
 252       use geometry_data, only: nres,rad2deg
 253       use energy_data
 254 !      implicit real*8 (a-h,o-z)
 255 !      include 'DIMENSIONS'
 256 !      include 'COMMON.IOUNITS'
 257 !      include 'COMMON.CSA'
 258 !      include 'COMMON.BANK'
 259 !      include 'COMMON.CHAIN'
 260 !      include 'COMMON.GEO'
 261 !      include 'COMMON.SBRIDGE'
 262 !el local variables
 263       integer :: nft,i,l,j,k,jlee
 264
 265       open(icsa_bank_reminimized,file=csa_bank_reminimized,&
 266         status="unknown")
 267       write (icsa_bank_reminimized,900) &
 268         jlee,nbank,nstep,nft,icycle,cutdif
 269       open (igeom,file=intname,status='UNKNOWN')
 270       do k=1,nbank
 271        write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
 272         bpncn(k),ibank(k)
 273        do j=1,numch
 274         do l=2,nres-1
 275          write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 276         enddo
 277        enddo
 278        if (nss.le.9) then
 279          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 280            nss,(ihpb(i),jhpb(i),i=1,nss)
 281        else
 282          write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
 283            nss,(ihpb(i),jhpb(i),i=1,9)
 284          write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
 285        endif
 286        write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
 287        write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
 288        write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
 289        write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
 290       enddo
 291       close(icsa_bank_reminimized)
 292       close(igeom)
 293
 294   200 format (8f10.4)
 295   850 format (10f8.3)
 296   900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
 297               i8,i10,i2,f15.5)
 298   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 299                ' %NC ',0pf5.2,i5)
 300
 301       return
 302       end subroutine write_bank_reminimized
 303 !-----------------------------------------------------------------------------
 304       subroutine read_bank(jlee,nft,cutdifr)
 305
 306       use csa_data
 307       use control_data, only: vdisulf
 308       use geometry_data, only: nres,deg2rad
 309       use energy_data
 310 !      implicit real*8 (a-h,o-z)
 311 !      include 'DIMENSIONS'
 312 !      include 'COMMON.IOUNITS'
 313 !      include 'COMMON.CSA'
 314 !      include 'COMMON.BANK'
 315 !      include 'COMMON.CHAIN'
 316 !      include 'COMMON.GEO'
 317 !      include 'COMMON.CONTROL'
 318 !      include 'COMMON.SBRIDGE'
 319       character(len=80) :: karta
 320 !      integer ilen
 321 !el      external ilen
 322 !el local variables
 323       integer :: nft,kk,k,l,i,j,jlee
 324       real(kind=8) :: cutdifr
 325
 326       open(icsa_bank,file=csa_bank,status="old")
 327        read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
 328        read (icsa_bank,902) nglob_csa, eglob_csa
 329 !      if(jleer.ne.jlee) then
 330 !        write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
 331 !    &   ' JLEE',jlee
 332 !        call mpi_abort(mpi_comm_world,ierror,ierrcode)
 333 !      endif
 334
 335       kk=0
 336       do k=1,nbank
 337         read (icsa_bank,'(a80)') karta
 338         write(iout,*) "READ_BANK: kk=",kk
 339         write(iout,*) karta
 340 !        if (index(karta,"*").gt.0) then
 341 !          write (iout,*) "***** Stars in bank ***** k=",k,
 342 !     &      " skipped"
 343 !          do j=1,numch
 344 !            do l=2,nres-1
 345 !              read (33,850) (rdummy,i=1,4)
 346 !            enddo
 347 !          enddo
 348 !        else
 349           kk=kk+1
 350           call reada(karta,"total E",bene(kk),1.0d20)
 351           call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
 352           call reada(karta,"%NC",bpncn(kk),0.0d0)
 353           read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
 354           goto 112
 355   111     ibank(kk)=0
 356   112     continue
 357           write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
 358             "%NC",bpncn(kk),ibank(kk)
 359 !          read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 360           if (vdisulf) then
 361             read (icsa_bank,'(101i4)') &
 362               bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 363             bvar_ns(kk)=ns-2*bvar_nss(kk)
 364             write(iout,*) 'read SSBOND',bvar_nss(kk),&
 365                           ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
 366 !d          write(iout,*) 'read CYS #free ', bvar_ns(kk)
 367             l=0
 368             do i=1,ns
 369              j=1
 370              do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
 371                        iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
 372                        j.le.bvar_nss(kk))
 373                 j=j+1
 374              enddo
 375              if (j.gt.bvar_nss(kk)) then
 376                l=l+1
 377                bvar_s(l,kk)=iss(i)
 378              endif
 379             enddo
 380 !d            write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
 381           endif
 382           do j=1,numch
 383             do l=2,nres-1
 384               read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
 385 !              write (iout,850) (bvar(i,l,j,kk),i=1,4)
 386               do i=1,4
 387                 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
 388               enddo ! l
 389             enddo ! l
 390           enddo ! j
 391 !        endif
 392       enddo ! k
 393
 394       if (kk.lt.nbank) nbank=kk
 395 !d      write (*,*) "read_bank ******************* kk",kk,
 396 !d     &  "nbank",nbank
 397 !d      do kk=1,nbank
 398 !d          print *,"kk=",kk
 399 !d          do j=1,numch
 400 !d            do l=2,nres-1
 401 !d              write (*,850) (bvar(i,l,j,kk),i=1,4)
 402 !d            enddo
 403 !d          enddo
 404 !d      enddo
 405
 406 !       do k=1,nbank
 407 !        read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
 408 !        do j=1,numch
 409 !         do l=2,nres-1
 410 !          read (33,850) (bvar(i,l,j,k),i=1,4)
 411 !          do i=1,4
 412 !           bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
 413 !          enddo
 414 !         enddo
 415 !        enddo
 416 !       enddo
 417       close(icsa_bank)
 418
 419   850 format (10f8.3)
 420   952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
 421   901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
 422   902 format (1x,11x,i4,12x,1pe14.5)
 423   953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
 424
 425       return
 426       end subroutine read_bank
 427 !-----------------------------------------------------------------------------
 428       subroutine write_bank1(jlee)
 429
 430       use csa_data
 431       use geometry_data, only: nres,rad2deg
 432 !      implicit real*8 (a-h,o-z)
 433 !      include 'DIMENSIONS'
 434 !      include 'COMMON.IOUNITS'
 435 !      include 'COMMON.CSA'
 436 !      include 'COMMON.BANK'
 437 !      include 'COMMON.CHAIN'
 438 !      include 'COMMON.GEO'
 439 !el local variables
 440       integer :: k,i,l,j,jlee
 441
 442 #if defined(AIX) || defined(PGI)
 443       open(icsa_bank1,file=csa_bank1,position="append")
 444 #else
 445       open(icsa_bank1,file=csa_bank1,access="append")
 446 #endif
 447       write (icsa_bank1,900) jlee,nbank,nstep,cutdif
 448       do k=1,nbank
 449        write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
 450        do j=1,numch
 451         do l=2,nres-1
 452          write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
 453         enddo
 454        enddo
 455       enddo
 456       close(icsa_bank1)
 457   850 format (10f8.3)
 458   900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
 459   952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
 460               ' %NC ',0pf5.2,i5)
 461
 462       return
 463       end subroutine write_bank1
 464 !-----------------------------------------------------------------------------
 465 ! cartprint.f
 466 !-----------------------------------------------------------------------------
 467 !      subroutine cartprint
 468
 469 !      use geometry_data, only: c
 470 !      use energy_data, only: itype
 471 !      implicit real*8 (a-h,o-z)
 472 !      include 'DIMENSIONS'
 473 !      include 'COMMON.CHAIN'
 474 !      include 'COMMON.INTERACT'
 475 !      include 'COMMON.NAMES'
 476 !      include 'COMMON.IOUNITS'
 477 !      integer :: i
 478
 479 !     write (iout,100)
 480 !      do i=1,nres
 481 !        write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
 482 !          c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
 483 !      enddo
 484 !  100 format (//'              alpha-carbon coordinates       ',&
 485 !                '     centroid coordinates'/ &
 486 !                '       ', 6X,'X',11X,'Y',11X,'Z',&
 487 !                                10X,'X',11X,'Y',11X,'Z')
 488 !  110 format (a,'(',i3,')',6f12.5)
 489 !      return
 490 !      end subroutine cartprint
 491 !-----------------------------------------------------------------------------
 492 ! dihed_cons.F
 493 !-----------------------------------------------------------------------------
 494       subroutine secstrp2dihc
 495
 496       use geometry_data
 497       use energy_data
 498 !      implicit real*8 (a-h,o-z)
 499 !      include 'DIMENSIONS'
 500 !      include 'COMMON.GEO'
 501 !      include 'COMMON.BOUNDS'
 502 !      include 'COMMON.CHAIN'
 503 !      include 'COMMON.TORCNSTR'
 504 !      include 'COMMON.IOUNITS'
 505 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 506 !el      COMMON/SECONDARYS/secstruc
 507       character(len=80) :: line
 508       logical :: errflag
 509 !el      external ilen
 510
 511 !el  local variables
 512       integer :: i,ii,lenpre
 513
 514       allocate(secstruc(nres))
 515
 516 !dr      call getenv_loc('SECPREDFIL',secpred)
 517       lenpre=ilen(prefix)
 518       secpred=prefix(:lenpre)//'.spred'
 519
 520 #if defined(WINIFL) || defined(WINPGI)
 521       open(isecpred,file=secpred,status='old',readonly,shared)
 522 #elif (defined CRAY) || (defined AIX)
 523       open(isecpred,file=secpred,status='old',action='read')
 524 #elif (defined G77)
 525       open(isecpred,file=secpred,status='old')
 526 #else
 527       open(isecpred,file=secpred,status='old',action='read')
 528 #endif
 529 ! read secondary structure prediction from JPRED here!
 530 !      read(isecpred,'(A80)',err=100,end=100) line
 531 !      read(line,'(f10.3)',err=110) ftors
 532        read(isecpred,'(f10.3)',err=110) ftors
 533
 534       write (iout,*) 'FTORS factor =',ftors
 535 ! initialize secstruc to any
 536        do i=1,nres
 537         secstruc(i) ='-'
 538       enddo
 539       ndih_constr=0
 540       ndih_nconstr=0
 541
 542       call read_secstr_pred(isecpred,iout,errflag)
 543       if (errflag) then
 544          write(iout,*)'There is a problem with the list of secondary-',&
 545            'structure prediction'
 546          goto 100
 547       endif
 548 ! 8/13/98 Set limits to generating the dihedral angles
 549       do i=1,nres
 550         phibound(1,i)=-pi
 551         phibound(2,i)=pi
 552       enddo
 553
 554       ii=0
 555       do i=1,nres
 556         if ( secstruc(i) .eq. 'H') then
 557 ! Helix restraints for this residue
 558            ii=ii+1
 559            idih_constr(ii)=i
 560            phi0(ii) = 45.0D0*deg2rad
 561            drange(ii)= 5.0D0*deg2rad
 562            phibound(1,i) = phi0(ii)-drange(ii)
 563            phibound(2,i) = phi0(ii)+drange(ii)
 564         else if (secstruc(i) .eq. 'E') then
 565 ! strand restraints for this residue
 566            ii=ii+1
 567            idih_constr(ii)=i
 568            phi0(ii) = 180.0D0*deg2rad
 569            drange(ii)= 5.0D0*deg2rad
 570            phibound(1,i) = phi0(ii)-drange(ii)
 571            phibound(2,i) = phi0(ii)+drange(ii)
 572         else
 573 ! no restraints for this residue
 574            ndih_nconstr=ndih_nconstr+1
 575            idih_nconstr(ndih_nconstr)=i
 576         endif
 577       enddo
 578       ndih_constr=ii
 579 !      deallocate(secstruc)
 580       return
 581 100   continue
 582       write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
 583 !      deallocate(secstruc)
 584       return
 585 110   continue
 586       write(iout,'(A20)')'Error reading FTORS'
 587 !      deallocate(secstruc)
 588       return
 589       end subroutine secstrp2dihc
 590 !-----------------------------------------------------------------------------
 591       subroutine read_secstr_pred(jin,jout,errors)
 592
 593 !      implicit real*8 (a-h,o-z)
 594 !      INCLUDE 'DIMENSIONS'
 595 !      include 'COMMON.IOUNITS'
 596 !      include 'COMMON.CHAIN'
 597 !el      character(len=1),dimension(nres) :: secstruc   !(maxres)
 598 !el      COMMON/SECONDARYS/secstruc
 599 !el      EXTERNAL ILEN
 600       character(len=80) :: line,line1   !,ucase
 601       logical :: errflag,errors,blankline
 602
 603 !el  local variables
 604       integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
 605             length_of_chain
 606       errors=.false.
 607       read (jin,'(a)') line
 608       write (jout,'(2a)') '> ',line(1:78)
 609       line1=ucase(line)
 610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
 611 ! correspond to the end-groups.  ADD to the secondary structure prediction "-" for the
 612 ! end-groups in the input file "*.spred"
 613
 614       iseq=1
 615       do while (index(line1,'$END').eq.0)
 616 ! Override commented lines.
 617          ipos=1
 618          blankline=.false.
 619          do while (.not.blankline)
 620             line1=' '
 621             call mykey(line,line1,ipos,blankline,errflag)
 622             if (errflag) write (jout,'(2a)') &
 623        'Error when reading sequence in line: ',line
 624             errors=errors .or. errflag
 625             if (.not. blankline .and. .not. errflag) then
 626                ipos1=2
 627                iend=ilen(line1)
 628 !el               if (iseq.le.maxres) then
 629                   if (line1(1:1).eq.'-' ) then
 630                      secstruc(iseq)=line1(1:1)
 631                   else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
 632                             ( ucase(line1(1:1)).eq.'H' ) ) then
 633                      secstruc(iseq)=ucase(line1(1:1))
 634                   else
 635                      errors=.true.
 636                      write (jout,1010) line1(1:1), iseq
 637                      goto 80
 638                   endif
 639 !el               else
 640 !el                  errors=.true.
 641 !el                  write (jout,1000) iseq,maxres
 642 !el                  goto 80
 643 !el               endif
 644                do while (ipos1.le.iend)
 645
 646                   iseq=iseq+1
 647                   il=1
 648                   ipos1=ipos1+1
 649 !el                  if (iseq.le.maxres) then
 650                      if (line1(ipos1-1:ipos1-1).eq.'-' ) then
 651                         secstruc(iseq)=line1(ipos1-1:ipos1-1)
 652                      else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
 653                            (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
 654                         secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
 655                      else
 656                         errors=.true.
 657                         write (jout,1010) line1(ipos1-1:ipos1-1), iseq
 658                         goto 80
 659                      endif
 660 !el                  else
 661 !el                     errors=.true.
 662 !el                     write (jout,1000) iseq,maxres
 663 !el                     goto 80
 664 !el                  endif
 665                enddo
 666                iseq=iseq+1
 667             endif
 668          enddo
 669          read (jin,'(a)') line
 670          write (jout,'(2a)') '> ',line(1:78)
 671          line1=ucase(line)
 672       enddo
 673
 674 !d    write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
 675
 676 !d check whether the found length of the chain is correct.
 677       length_of_chain=iseq-1
 678       if (length_of_chain .ne. nres) then
 679 !        errors=.true.
 680         write (jout,'(a,i4,a,i4,a)') &
 681        'Error: the number of labels specified in $SEC_STRUC_PRED (' &
 682        ,length_of_chain,') does not match with the number of residues (' &
 683        ,nres,').'
 684       endif
 685    80 continue
 686
 687  1000 format('Error - the number of residues (',i4,&
 688        ') has exceeded maximum (',i4,').')
 689  1010 format ('Error - unrecognized secondary structure label',a4,&
 690        ' in position',i4)
 691       return
 692       end subroutine read_secstr_pred
 693 !#endif
 694 !-----------------------------------------------------------------------------
 695 ! parmread.F
 696 !-----------------------------------------------------------------------------
 697       subroutine parmread
 698
 699       use geometry_data
 700       use energy_data
 701       use control_data, only:maxterm !,maxtor
 702       use MD_data
 703       use MPI_data
 704 !el      use map_data
 705       use control, only: getenv_loc
 706 !
 707 ! Read the parameters of the probability distributions of the virtual-bond
 708 ! valence angles and the side chains and energy parameters.
 709 !
 710 ! Important! Energy-term weights ARE NOT read here; they are read from the
 711 ! main input file instead, because NO defaults have yet been set for these
 712 ! parameters.
 713 !
 714 !      implicit real*8 (a-h,o-z)
 715 !      include 'DIMENSIONS'
 716 #ifdef MPI
 717       include "mpif.h"
 718       integer :: IERROR
 719 #endif
 720 !      include 'COMMON.IOUNITS'
 721 !      include 'COMMON.CHAIN'
 722 !      include 'COMMON.INTERACT'
 723 !      include 'COMMON.GEO'
 724 !      include 'COMMON.LOCAL'
 725 !      include 'COMMON.TORSION'
 726 !      include 'COMMON.SCCOR'
 727 !      include 'COMMON.SCROT'
 728 !      include 'COMMON.FFIELD'
 729 !      include 'COMMON.NAMES'
 730 !      include 'COMMON.SBRIDGE'
 731 !      include 'COMMON.MD'
 732 !      include 'COMMON.SETUP'
 733       character(len=1) :: t1,t2,t3
 734       character(len=1) :: onelett(4) = (/"G","A","P","D"/)
 735       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
 736       logical :: lprint,LaTeX
 737       real(kind=8),dimension(3,3,maxlob) :: blower      !(3,3,maxlob)
 738       real(kind=8),dimension(13) :: b
 739       character(len=3) :: lancuch       !,ucase
 740 !el  local variables
 741       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
 742       integer :: maxinter,junk,kk,ii,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         write (iout,'(/a/)') 'Torsional constants:'
1901         do i=1,nsccortyp
1902           do j=1,nsccortyp
1903             write (iout,*) 'ityp',i,' jtyp',j
1904             write (iout,*) 'Fourier constants'
1905             do k=1,nterm_sccor(i,j)
1906       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1907             enddo
1908             write (iout,*) 'Lorenz constants'
1909             do k=1,nlor_sccor(i,j)
1910               write (iout,'(3(1pe15.5))') &
1911                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1912             enddo
1913           enddo
1914         enddo
1915       endif
1916
1917 !
1918 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1919 !         interaction energy of the Gly, Ala, and Pro prototypes.
1920 !
1921
1922       if (lprint) then
1923         write (iout,*)
1924         write (iout,*) "Coefficients of the cumulants"
1925       endif
1926       read (ifourier,*) nloctyp
1927 !write(iout,*) "nloctyp",nloctyp
1928 !el from module energy-------
1929       allocate(b1(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1930       allocate(b2(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1931       allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1932       allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1933       allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1934       allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1935       allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1936       allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1937 ! el
1938         b1(1,:)=0.0d0
1939         b1(2,:)=0.0d0
1940 !--------------------------------
1941
1942       do i=0,nloctyp-1
1943         read (ifourier,*,end=115,err=115)
1944         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1945         if (lprint) then
1946           write (iout,*) 'Type',i
1947           write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1948         endif
1949         B1(1,i)  = b(3)
1950         B1(2,i)  = b(5)
1951         B1(1,-i) = b(3)
1952         B1(2,-i) = -b(5)
1953 !        b1(1,i)=0.0d0
1954 !        b1(2,i)=0.0d0
1955         B1tilde(1,i) = b(3)
1956         B1tilde(2,i) =-b(5)
1957         B1tilde(1,-i) =-b(3)
1958         B1tilde(2,-i) =b(5)
1959 !        b1tilde(1,i)=0.0d0
1960 !        b1tilde(2,i)=0.0d0
1961         B2(1,i)  = b(2)
1962         B2(2,i)  = b(4)
1963         B2(1,-i)  =b(2)
1964         B2(2,-i)  =-b(4)
1965
1966 !        b2(1,i)=0.0d0
1967 !        b2(2,i)=0.0d0
1968         CC(1,1,i)= b(7)
1969         CC(2,2,i)=-b(7)
1970         CC(2,1,i)= b(9)
1971         CC(1,2,i)= b(9)
1972         CC(1,1,-i)= b(7)
1973         CC(2,2,-i)=-b(7)
1974         CC(2,1,-i)=-b(9)
1975         CC(1,2,-i)=-b(9)
1976 !        CC(1,1,i)=0.0d0
1977 !        CC(2,2,i)=0.0d0
1978 !        CC(2,1,i)=0.0d0
1979 !        CC(1,2,i)=0.0d0
1980         Ctilde(1,1,i)=b(7)
1981         Ctilde(1,2,i)=b(9)
1982         Ctilde(2,1,i)=-b(9)
1983         Ctilde(2,2,i)=b(7)
1984         Ctilde(1,1,-i)=b(7)
1985         Ctilde(1,2,-i)=-b(9)
1986         Ctilde(2,1,-i)=b(9)
1987         Ctilde(2,2,-i)=b(7)
1988
1989 !        Ctilde(1,1,i)=0.0d0
1990 !        Ctilde(1,2,i)=0.0d0
1991 !        Ctilde(2,1,i)=0.0d0
1992 !        Ctilde(2,2,i)=0.0d0
1993         DD(1,1,i)= b(6)
1994         DD(2,2,i)=-b(6)
1995         DD(2,1,i)= b(8)
1996         DD(1,2,i)= b(8)
1997         DD(1,1,-i)= b(6)
1998         DD(2,2,-i)=-b(6)
1999         DD(2,1,-i)=-b(8)
2000         DD(1,2,-i)=-b(8)
2001 !        DD(1,1,i)=0.0d0
2002 !        DD(2,2,i)=0.0d0
2003 !        DD(2,1,i)=0.0d0
2004 !        DD(1,2,i)=0.0d0
2005         Dtilde(1,1,i)=b(6)
2006         Dtilde(1,2,i)=b(8)
2007         Dtilde(2,1,i)=-b(8)
2008         Dtilde(2,2,i)=b(6)
2009         Dtilde(1,1,-i)=b(6)
2010         Dtilde(1,2,-i)=-b(8)
2011         Dtilde(2,1,-i)=b(8)
2012         Dtilde(2,2,-i)=b(6)
2013
2014 !        Dtilde(1,1,i)=0.0d0
2015 !        Dtilde(1,2,i)=0.0d0
2016 !        Dtilde(2,1,i)=0.0d0
2017 !        Dtilde(2,2,i)=0.0d0
2018         EE(1,1,i)= b(10)+b(11)
2019         EE(2,2,i)=-b(10)+b(11)
2020         EE(2,1,i)= b(12)-b(13)
2021         EE(1,2,i)= b(12)+b(13)
2022         EE(1,1,-i)= b(10)+b(11)
2023         EE(2,2,-i)=-b(10)+b(11)
2024         EE(2,1,-i)=-b(12)+b(13)
2025         EE(1,2,-i)=-b(12)-b(13)
2026
2027 !        ee(1,1,i)=1.0d0
2028 !        ee(2,2,i)=1.0d0
2029 !        ee(2,1,i)=0.0d0
2030 !        ee(1,2,i)=0.0d0
2031 !        ee(2,1,i)=ee(1,2,i)
2032       enddo
2033       if (lprint) then
2034       do i=1,nloctyp
2035         write (iout,*) 'Type',i
2036         write (iout,*) 'B1'
2037         write(iout,*) B1(1,i),B1(2,i)
2038         write (iout,*) 'B2'
2039         write(iout,*) B2(1,i),B2(2,i)
2040         write (iout,*) 'CC'
2041         do j=1,2
2042           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2043         enddo
2044         write(iout,*) 'DD'
2045         do j=1,2
2046           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2047         enddo
2048         write(iout,*) 'EE'
2049         do j=1,2
2050           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2051         enddo
2052       enddo
2053       endif
2054 !
2055 ! Read electrostatic-interaction parameters
2056 !
2057
2058       if (lprint) then
2059         write (iout,*)
2060         write (iout,'(/a)') 'Electrostatic interaction constants:'
2061         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2062                   'IT','JT','APP','BPP','AEL6','AEL3'
2063       endif
2064       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2065       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2066       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2067       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2068       close (ielep)
2069       do i=1,2
2070         do j=1,2
2071         rri=rpp(i,j)**6
2072         app (i,j)=epp(i,j)*rri*rri
2073         bpp (i,j)=-2.0D0*epp(i,j)*rri
2074         ael6(i,j)=elpp6(i,j)*4.2D0**6
2075         ael3(i,j)=elpp3(i,j)*4.2D0**3
2076 !        lprint=.true.
2077         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2078                           ael6(i,j),ael3(i,j)
2079 !        lprint=.false.
2080         enddo
2081       enddo
2082 !
2083 ! Read side-chain interaction parameters.
2084 !
2085 !el from module energy - COMMON.INTERACT-------
2086       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2087       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2088       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2089
2090       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2091       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2092       allocate(epslip(ntyp,ntyp))
2093       augm(:,:)=0.0D0
2094       chip(:)=0.0D0
2095       alp(:)=0.0D0
2096       sigma0(:)=0.0D0
2097       sigii(:)=0.0D0
2098       rr0(:)=0.0D0
2099
2100 !--------------------------------
2101
2102       read (isidep,*,end=117,err=117) ipot,expon
2103       if (ipot.lt.1 .or. ipot.gt.5) then
2104         write (iout,'(2a)') 'Error while reading SC interaction',&
2105                      'potential file - unknown potential type.'
2106 #ifdef MPI
2107         call MPI_Finalize(Ierror)
2108 #endif
2109         stop
2110       endif
2111       expon2=expon/2
2112       if(me.eq.king) &
2113        write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2114        ', exponents are ',expon,2*expon
2115 !      goto (10,20,30,30,40) ipot
2116       select case(ipot)
2117 !----------------------- LJ potential ---------------------------------
2118        case (1)
2119 !   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2120          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2121            (sigma0(i),i=1,ntyp)
2122         if (lprint) then
2123           write (iout,'(/a/)') 'Parameters of the LJ potential:'
2124           write (iout,'(a/)') 'The epsilon array:'
2125           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2126           write (iout,'(/a)') 'One-body parameters:'
2127           write (iout,'(a,4x,a)') 'residue','sigma'
2128           write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2129         endif
2130 !      goto 50
2131 !----------------------- LJK potential --------------------------------
2132        case(2)
2133 !   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2134          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2135           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2136         if (lprint) then
2137           write (iout,'(/a/)') 'Parameters of the LJK potential:'
2138           write (iout,'(a/)') 'The epsilon array:'
2139           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2140           write (iout,'(/a)') 'One-body parameters:'
2141           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
2142           write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2143                 i=1,ntyp)
2144         endif
2145 !      goto 50
2146 !---------------------- GB or BP potential -----------------------------
2147        case(3:4)
2148 !   30 do i=1,ntyp
2149         do i=1,ntyp
2150          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2151         enddo
2152         read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2153         read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2154         read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2155         read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2156         do i=1,ntyp
2157          read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2158         enddo
2159
2160 ! For the GB potential convert sigma'**2 into chi'
2161         if (ipot.eq.4) then
2162           do i=1,ntyp
2163             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2164           enddo
2165         endif
2166         if (lprint) then
2167           write (iout,'(/a/)') 'Parameters of the BP potential:'
2168           write (iout,'(a/)') 'The epsilon array:'
2169           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2170           write (iout,'(/a)') 'One-body parameters:'
2171           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
2172                '    chip  ','    alph  '
2173           write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2174                              chip(i),alp(i),i=1,ntyp)
2175         endif
2176 !      goto 50
2177 !--------------------- GBV potential -----------------------------------
2178        case(5)
2179 !   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2180         read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2181           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2182           (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2183         if (lprint) then
2184           write (iout,'(/a/)') 'Parameters of the GBV potential:'
2185           write (iout,'(a/)') 'The epsilon array:'
2186           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2187           write (iout,'(/a)') 'One-body parameters:'
2188           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
2189               's||/s_|_^2','    chip  ','    alph  '
2190           write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2191                    sigii(i),chip(i),alp(i),i=1,ntyp)
2192         endif
2193        case default
2194         write(iout,*)"Wrong ipot"
2195 !   50 continue
2196       end select
2197       continue
2198       close (isidep)
2199
2200 !-----------------------------------------------------------------------
2201 ! Calculate the "working" parameters of SC interactions.
2202
2203 !el from module energy - COMMON.INTERACT-------
2204       allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2205       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2206       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2207       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2208         dcavtub(ntyp1))
2209       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2210         tubetranene(ntyp1))
2211       aa_aq(:,:)=0.0D0
2212       bb_aq(:,:)=0.0D0
2213       aa_lip(:,:)=0.0D0
2214       bb_lip(:,:)=0.0D0
2215       chi(:,:)=0.0D0
2216       sigma(:,:)=0.0D0
2217       r0(:,:)=0.0D0
2218       acavtub(:)=0.0d0
2219       bcavtub(:)=0.0d0
2220       ccavtub(:)=0.0d0
2221       dcavtub(:)=0.0d0
2222       sc_aa_tube_par(:)=0.0d0
2223       sc_bb_tube_par(:)=0.0d0
2224
2225 !--------------------------------
2226
2227       do i=2,ntyp
2228         do j=1,i-1
2229           eps(i,j)=eps(j,i)
2230           epslip(i,j)=epslip(j,i)
2231         enddo
2232       enddo
2233       do i=1,ntyp
2234         do j=i,ntyp
2235           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2236           sigma(j,i)=sigma(i,j)
2237           rs0(i,j)=dwa16*sigma(i,j)
2238           rs0(j,i)=rs0(i,j)
2239         enddo
2240       enddo
2241       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2242        'Working parameters of the SC interactions:',&
2243        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
2244        '  chi1   ','   chi2   '
2245       do i=1,ntyp
2246         do j=i,ntyp
2247           epsij=eps(i,j)
2248           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2249             rrij=sigma(i,j)
2250           else
2251             rrij=rr0(i)+rr0(j)
2252           endif
2253           r0(i,j)=rrij
2254           r0(j,i)=rrij
2255           rrij=rrij**expon
2256           epsij=eps(i,j)
2257           sigeps=dsign(1.0D0,epsij)
2258           epsij=dabs(epsij)
2259           aa_aq(i,j)=epsij*rrij*rrij
2260           bb_aq(i,j)=-sigeps*epsij*rrij
2261           aa_aq(j,i)=aa_aq(i,j)
2262           bb_aq(j,i)=bb_aq(i,j)
2263           epsijlip=epslip(i,j)
2264           sigeps=dsign(1.0D0,epsijlip)
2265           epsijlip=dabs(epsijlip)
2266           aa_lip(i,j)=epsijlip*rrij*rrij
2267           bb_lip(i,j)=-sigeps*epsijlip*rrij
2268           aa_lip(j,i)=aa_lip(i,j)
2269           bb_lip(j,i)=bb_lip(i,j)
2270 !C          write(iout,*) aa_lip
2271           if (ipot.gt.2) then
2272             sigt1sq=sigma0(i)**2
2273             sigt2sq=sigma0(j)**2
2274             sigii1=sigii(i)
2275             sigii2=sigii(j)
2276             ratsig1=sigt2sq/sigt1sq
2277             ratsig2=1.0D0/ratsig1
2278             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2279             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2280             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2281           else
2282             rsum_max=sigma(i,j)
2283           endif
2284 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2285             sigmaii(i,j)=rsum_max
2286             sigmaii(j,i)=rsum_max
2287 !         else
2288 !           sigmaii(i,j)=r0(i,j)
2289 !           sigmaii(j,i)=r0(i,j)
2290 !         endif
2291 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2292           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2293             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2294             augm(i,j)=epsij*r_augm**(2*expon)
2295 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2296             augm(j,i)=augm(i,j)
2297           else
2298             augm(i,j)=0.0D0
2299             augm(j,i)=0.0D0
2300           endif
2301           if (lprint) then
2302             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2303             restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2304             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2305           endif
2306         enddo
2307       enddo
2308
2309       allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2310       allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2311       allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2312       allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2313       allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2314       allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2315       allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2316       allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2317       allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2318       allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2319       allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2320       allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2321       allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2322       allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2323       allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2324       allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2325
2326 !      augm(:,:)=0.0D0
2327 !      chip(:)=0.0D0
2328 !      alp(:)=0.0D0
2329 !      sigma0(:)=0.0D0
2330 !      sigii(:)=0.0D0
2331 !      rr0(:)=0.0D0
2332
2333       read (isidep_nucl,*) ipot_nucl
2334 !      print *,"TU?!",ipot_nucl
2335       if (ipot_nucl.eq.1) then
2336         do i=1,ntyp_molec(2)
2337           do j=i,ntyp_molec(2)
2338             read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2339             elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2340           enddo
2341         enddo
2342       else
2343         do i=1,ntyp_molec(2)
2344           do j=i,ntyp_molec(2)
2345             read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2346                chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2347                elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2348           enddo
2349         enddo
2350       endif
2351 !      rpp(1,1)=2**(1.0/6.0)*5.16158
2352       do i=1,ntyp_molec(2)
2353         do j=i,ntyp_molec(2)
2354           rrij=sigma_nucl(i,j)
2355           r0_nucl(i,j)=rrij
2356           r0_nucl(j,i)=rrij
2357           rrij=rrij**expon
2358           epsij=4*eps_nucl(i,j)
2359           sigeps=dsign(1.0D0,epsij)
2360           epsij=dabs(epsij)
2361           aa_nucl(i,j)=epsij*rrij*rrij
2362           bb_nucl(i,j)=-sigeps*epsij*rrij
2363           ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2364           ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2365           ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2366           ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2367           sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2368          2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2369         enddo
2370         do j=1,i-1
2371           aa_nucl(i,j)=aa_nucl(j,i)
2372           bb_nucl(i,j)=bb_nucl(j,i)
2373           ael3_nucl(i,j)=ael3_nucl(j,i)
2374           ael6_nucl(i,j)=ael6_nucl(j,i)
2375           ael63_nucl(i,j)=ael63_nucl(j,i)
2376           ael32_nucl(i,j)=ael32_nucl(j,i)
2377           elpp3_nucl(i,j)=elpp3_nucl(j,i)
2378           elpp6_nucl(i,j)=elpp6_nucl(j,i)
2379           elpp63_nucl(i,j)=elpp63_nucl(j,i)
2380           elpp32_nucl(i,j)=elpp32_nucl(j,i)
2381           eps_nucl(i,j)=eps_nucl(j,i)
2382           sigma_nucl(i,j)=sigma_nucl(j,i)
2383           sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2384         enddo
2385       enddo
2386
2387       write(iout,*) "tube param"
2388       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2389       ccavtubpep,dcavtubpep,tubetranenepep
2390       sigmapeptube=sigmapeptube**6
2391       sigeps=dsign(1.0D0,epspeptube)
2392       epspeptube=dabs(epspeptube)
2393       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2394       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2395       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2396       do i=1,ntyp
2397        read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2398       ccavtub(i),dcavtub(i),tubetranene(i)
2399        sigmasctube=sigmasctube**6
2400        sigeps=dsign(1.0D0,epssctube)
2401        epssctube=dabs(epssctube)
2402        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2403        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2404       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2405       enddo
2406 !-----------------READING SC BASE POTENTIALS-----------------------------
2407       allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2408       allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2409       allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2410       allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2411       allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2412       allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2413       allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2414       allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2415       allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2416       allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2417       allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2418       allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2419       allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2420       allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2421       allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2422       allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2423
2424
2425       do i=1,ntyp_molec(1)
2426        do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2427         write (*,*) "Im in ", i, " ", j
2428         read(isidep_scbase,*) &
2429         eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2430         chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2431          write(*,*) "eps",eps_scbase(i,j)
2432         read(isidep_scbase,*) &
2433        (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2434        chis_scbase(i,j,1),chis_scbase(i,j,2)
2435         read(isidep_scbase,*) &
2436        dhead_scbasei(i,j), &
2437        dhead_scbasej(i,j), &
2438        rborn_scbasei(i,j),rborn_scbasej(j,i)
2439         read(isidep_scbase,*) &
2440        (wdipdip_scbase(k,i,j),k=1,3), &
2441        (wqdip_scbase(k,i,j),k=1,2)
2442         read(isidep_scbase,*) &
2443        alphapol_scbase(i,j), &
2444        epsintab_scbase(i,j)
2445        END DO
2446       END DO
2447       allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2448       allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2449
2450       do i=1,ntyp_molec(1)
2451        do j=1,ntyp_molec(2)-1
2452           epsij=eps_scbase(i,j)
2453           rrij=sigma(i,j)
2454 !          r0(i,j)=rrij
2455 !          r0(j,i)=rrij
2456           rrij=rrij**expon
2457 !          epsij=eps(i,j)
2458           sigeps=dsign(1.0D0,epsij)
2459           epsij=dabs(epsij)
2460           aa_scbase(i,j)=epsij*rrij*rrij
2461           bb_scbase(i,j)=-sigeps*epsij*rrij
2462         enddo
2463        enddo
2464
2465       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2466       bad(:,:)=0.0D0
2467
2468 #ifdef OLDSCP
2469 !
2470 ! Define the SC-p interaction constants (hard-coded; old style)
2471 !
2472       do i=1,ntyp
2473 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2474 ! helix formation)
2475 !       aad(i,1)=0.3D0*4.0D0**12
2476 ! Following line for constants currently implemented
2477 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2478         aad(i,1)=1.5D0*4.0D0**12
2479 !       aad(i,1)=0.17D0*5.6D0**12
2480         aad(i,2)=aad(i,1)
2481 ! "Soft" SC-p repulsion
2482         bad(i,1)=0.0D0
2483 ! Following line for constants currently implemented
2484 !       aad(i,1)=0.3D0*4.0D0**6
2485 ! "Hard" SC-p repulsion
2486         bad(i,1)=3.0D0*4.0D0**6
2487 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2488         bad(i,2)=bad(i,1)
2489 !       aad(i,1)=0.0D0
2490 !       aad(i,2)=0.0D0
2491 !       bad(i,1)=1228.8D0
2492 !       bad(i,2)=1228.8D0
2493       enddo
2494 #else
2495 !
2496 ! 8/9/01 Read the SC-p interaction constants from file
2497 !
2498       do i=1,ntyp
2499         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2500       enddo
2501       do i=1,ntyp
2502         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2503         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2504         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2505         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2506       enddo
2507 !      lprint=.true.
2508       if (lprint) then
2509         write (iout,*) "Parameters of SC-p interactions:"
2510         do i=1,ntyp
2511           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2512            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2513         enddo
2514       endif
2515 !      lprint=.false.
2516 #endif
2517       allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2518
2519       do i=1,ntyp_molec(2)
2520         read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2521       enddo
2522       do i=1,ntyp_molec(2)
2523         aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2524         bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2525       enddo
2526       r0pp=1.12246204830937298142*5.16158
2527       epspp=4.95713/4
2528       AEES=108.661
2529       BEES=0.433246
2530
2531 !
2532 ! Define the constants of the disulfide bridge
2533 !
2534       ebr=-5.50D0
2535 !
2536 ! Old arbitrary potential - commented out.
2537 !
2538 !      dbr= 4.20D0
2539 !      fbr= 3.30D0
2540 !
2541 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2542 ! energy surface of diethyl disulfide.
2543 ! A. Liwo and U. Kozlowska, 11/24/03
2544 !
2545       D0CM = 3.78d0
2546       AKCM = 15.1d0
2547       AKTH = 11.0d0
2548       AKCT = 12.0d0
2549       V1SS =-1.08d0
2550       V2SS = 7.61d0
2551       V3SS = 13.7d0
2552 !      akcm=0.0d0
2553 !      akth=0.0d0
2554 !      akct=0.0d0
2555 !      v1ss=0.0d0
2556 !      v2ss=0.0d0
2557 !      v3ss=0.0d0
2558
2559       if(me.eq.king) then
2560       write (iout,'(/a)') "Disulfide bridge parameters:"
2561       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2562       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2563       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2564       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2565         ' v3ss:',v3ss
2566       endif
2567       return
2568   111 write (iout,*) "Error reading bending energy parameters."
2569       goto 999
2570   112 write (iout,*) "Error reading rotamer energy parameters."
2571       goto 999
2572   113 write (iout,*) "Error reading torsional energy parameters."
2573       goto 999
2574   114 write (iout,*) "Error reading double torsional energy parameters."
2575       goto 999
2576   115 write (iout,*) &
2577         "Error reading cumulant (multibody energy) parameters."
2578       goto 999
2579   116 write (iout,*) "Error reading electrostatic energy parameters."
2580       goto 999
2581   117 write (iout,*) "Error reading side chain interaction parameters."
2582       goto 999
2583   118 write (iout,*) "Error reading SCp interaction parameters."
2584       goto 999
2585   119 write (iout,*) "Error reading SCCOR parameters"
2586   999 continue
2587 #ifdef MPI
2588       call MPI_Finalize(Ierror)
2589 #endif
2590       stop
2591       return
2592       end subroutine parmread
2593 #endif
2594 !-----------------------------------------------------------------------------
2595 ! printmat.f
2596 !-----------------------------------------------------------------------------
2597       subroutine printmat(ldim,m,n,iout,key,a)
2598
2599       integer :: n,ldim
2600       character(len=3),dimension(n) :: key
2601       real(kind=8),dimension(ldim,n) :: a
2602 !el local variables
2603       integer :: i,j,k,m,iout,nlim
2604
2605       do 1 i=1,n,8
2606       nlim=min0(i+7,n)
2607       write (iout,1000) (key(k),k=i,nlim)
2608       write (iout,1020)
2609  1000 format (/5x,8(6x,a3))
2610  1020 format (/80(1h-)/)
2611       do 2 j=1,n
2612       write (iout,1010) key(j),(a(j,k),k=i,nlim)
2613     2 continue
2614     1 continue
2615  1010 format (a3,2x,8(f9.4))
2616       return
2617       end subroutine printmat
2618 !-----------------------------------------------------------------------------
2619 ! readpdb.F
2620 !-----------------------------------------------------------------------------
2621       subroutine readpdb
2622 ! Read the PDB file and convert the peptide geometry into virtual-chain
2623 ! geometry.
2624       use geometry_data
2625       use energy_data, only: itype,istype
2626       use control_data
2627       use compare_data
2628       use MPI_data
2629       use control, only: rescode,sugarcode
2630 !      implicit real*8 (a-h,o-z)
2631 !      include 'DIMENSIONS'
2632 !      include 'COMMON.LOCAL'
2633 !      include 'COMMON.VAR'
2634 !      include 'COMMON.CHAIN'
2635 !      include 'COMMON.INTERACT'
2636 !      include 'COMMON.IOUNITS'
2637 !      include 'COMMON.GEO'
2638 !      include 'COMMON.NAMES'
2639 !      include 'COMMON.CONTROL'
2640 !      include 'COMMON.DISTFIT'
2641 !      include 'COMMON.SETUP'
2642       integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2643 !        ishift_pdb
2644       logical :: lprn=.true.,fail
2645       real(kind=8),dimension(3) :: e1,e2,e3
2646       real(kind=8) :: dcj,efree_temp
2647       character(len=3) :: seq,res
2648       character(len=5) :: atom
2649       character(len=80) :: card
2650       real(kind=8),dimension(3,20) :: sccor
2651       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin       !rescode,
2652       integer :: isugar,molecprev,firstion
2653       character*1 :: sugar
2654       real(kind=8) :: cou
2655       real(kind=8),dimension(3) :: ccc
2656 !el local varables
2657       integer,dimension(2,maxres/3) :: hfrag_alloc
2658       integer,dimension(4,maxres/3) :: bfrag_alloc
2659       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2660       real(kind=8),dimension(:,:), allocatable  :: c_temporary
2661       integer,dimension(:,:) , allocatable  :: itype_temporary
2662       integer,dimension(:),allocatable :: istype_temp
2663       efree_temp=0.0d0
2664       ibeg=1
2665       ishift1=0
2666       ishift=0
2667       molecule=1
2668       counter=0
2669 !      write (2,*) "UNRES_PDB",unres_pdb
2670       ires=0
2671       ires_old=0
2672       nres=0
2673       iii=0
2674       lsecondary=.false.
2675       nhfrag=0
2676       nbfrag=0
2677 !-----------------------------
2678       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2679       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2680
2681       do i=1,100000
2682         read (ipdbin,'(a80)',end=10) card
2683 !       write (iout,'(a)') card
2684         if (card(:5).eq.'HELIX') then
2685           nhfrag=nhfrag+1
2686           lsecondary=.true.
2687           read(card(22:25),*) hfrag(1,nhfrag)
2688           read(card(34:37),*) hfrag(2,nhfrag)
2689         endif
2690         if (card(:5).eq.'SHEET') then
2691           nbfrag=nbfrag+1
2692           lsecondary=.true.
2693           read(card(24:26),*) bfrag(1,nbfrag)
2694           read(card(35:37),*) bfrag(2,nbfrag)
2695 !rc----------------------------------------
2696 !rc  to be corrected !!!
2697           bfrag(3,nbfrag)=bfrag(1,nbfrag)
2698           bfrag(4,nbfrag)=bfrag(2,nbfrag)
2699 !rc----------------------------------------
2700         endif
2701         if (card(:3).eq.'END') then
2702           goto 10
2703         else if (card(:3).eq.'TER') then
2704 ! End current chain
2705           ires_old=ires+2
2706           ishift=ishift+1
2707           ishift1=ishift1+1
2708           itype(ires_old,molecule)=ntyp1_molec(molecule)
2709           itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2710           nres_molec(molecule)=nres_molec(molecule)+2
2711           ibeg=2
2712 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2713           if (unres_pdb) then
2714             do j=1,3
2715               dc(j,ires)=sccor(j,iii)
2716             enddo
2717           else
2718             call sccenter(ires,iii,sccor)
2719 !          iii=0
2720           endif
2721           iii=0
2722         endif
2723 ! Read free energy
2724         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2725 ! Fish out the ATOM cards.
2726 !        write(iout,*) 'card',card(1:20)
2727         if (index(card(1:4),'ATOM').gt.0) then
2728           read (card(12:16),*) atom
2729 !          write (iout,*) "! ",atom," !",ires
2730 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2731           read (card(23:26),*) ires
2732           read (card(18:20),'(a3)') res
2733 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2734 !     &      " ires_old",ires_old
2735 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2736 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2737           if (ires-ishift+ishift1.ne.ires_old) then
2738 ! Calculate the CM of the preceding residue.
2739 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2740             if (ibeg.eq.0) then
2741 !              write (iout,*) "Calculating sidechain center iii",iii
2742               if (unres_pdb) then
2743                 do j=1,3
2744                   dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2745                 enddo
2746               else
2747                 call sccenter(ires_old,iii,sccor)
2748               endif
2749               iii=0
2750             endif
2751 ! Start new residue.
2752             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2753               ires=ires_old
2754               cycle
2755             else if (ibeg.eq.1) then
2756               write (iout,*) "BEG ires",ires
2757               ishift=ires-1
2758               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2759                 ishift=ishift-1
2760                 itype(1,1)=ntyp1
2761                 nres_molec(molecule)=nres_molec(molecule)+1
2762               endif
2763               ires=ires-ishift+ishift1
2764               ires_old=ires
2765 !              write (iout,*) "ishift",ishift," ires",ires,&
2766 !               " ires_old",ires_old
2767               ibeg=0
2768             else if (ibeg.eq.2) then
2769 ! Start a new chain
2770               ishift=-ires_old+ires-1 !!!!!
2771               ishift1=ishift1-1    !!!!!
2772 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2773               ires=ires-ishift+ishift1
2774 !              print *,ires,ishift,ishift1
2775               ires_old=ires
2776               ibeg=0
2777             else
2778               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2779               ires=ires-ishift+ishift1
2780               ires_old=ires
2781             endif
2782 !            print *,'atom',ires,atom
2783             if (res.eq.'ACE' .or. res.eq.'NHE') then
2784               itype(ires,1)=10
2785             else
2786              if (atom.eq.'CA  '.or.atom.eq.'N   ') then
2787              molecule=1
2788               itype(ires,molecule)=rescode(ires,res,0,molecule)
2789               firstion=0
2790 !              nres_molec(molecule)=nres_molec(molecule)+1
2791             else
2792              molecule=2
2793               itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2794 !              nres_molec(molecule)=nres_molec(molecule)+1
2795              read (card(19:19),'(a1)') sugar
2796              isugar=sugarcode(sugar,ires)
2797 !            if (ibeg.eq.1) then
2798 !              istype(1)=isugar
2799 !            else
2800               istype(ires)=isugar
2801 !              print *,"ires=",ires,istype(ires)
2802 !            endif
2803
2804             endif
2805             endif
2806           else
2807             ires=ires-ishift+ishift1
2808           endif
2809 !          write (iout,*) "ires_old",ires_old," ires",ires
2810           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2811 !            ishift1=ishift1+1
2812           endif
2813 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2814           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2815              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2816             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2817 !              print *,ires,ishift,ishift1
2818 !            write (iout,*) "backbone ",atom
2819 #ifdef DEBUG
2820             write (iout,'(2i3,2x,a,3f8.3)') &
2821             ires,itype(ires,1),res,(c(j,ires),j=1,3)
2822 #endif
2823             iii=iii+1
2824               nres_molec(molecule)=nres_molec(molecule)+1
2825             do j=1,3
2826               sccor(j,iii)=c(j,ires)
2827             enddo
2828           else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2829                atom.eq."C2'" .or. atom.eq."C3'" &
2830                .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2831             read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2832 !c            write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2833 !              print *,ires,ishift,ishift1
2834             counter=counter+1
2835 !            iii=iii+1
2836 !            do j=1,3
2837 !              sccor(j,iii)=c(j,ires)
2838 !            enddo
2839             do j=1,3
2840               c(j,ires)=c(j,ires)+ccc(j)/5.0
2841             enddo
2842              print *,counter,molecule
2843              if (counter.eq.5) then
2844 !            iii=iii+1
2845               nres_molec(molecule)=nres_molec(molecule)+1
2846               firstion=0
2847 !            do j=1,3
2848 !              sccor(j,iii)=c(j,ires)
2849 !            enddo
2850              counter=0
2851            endif
2852 !            print *, "ATOM",atom(1:3)
2853           else if (atom.eq."C5'") then
2854              read (card(19:19),'(a1)') sugar
2855              isugar=sugarcode(sugar,ires)
2856             if (ibeg.eq.1) then
2857               istype(1)=isugar
2858             else
2859               istype(ires)=isugar
2860 !              print *,ires,istype(ires)
2861             endif
2862             if (unres_pdb) then
2863               read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2864             else
2865               iii=iii+1
2866               read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2867             endif
2868 !            write (*,*) card(23:27),ires,itype(ires,1)
2869           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2870                    atom.ne.'N' .and. atom.ne.'C' .and. &
2871                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2872                    atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2873                    .and. atom.ne.'P  '.and. &
2874                   atom(1:1).ne.'H' .and. &
2875                   atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2876                   ) then
2877 !            write (iout,*) "sidechain ",atom
2878 !            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2879                  if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2880 !                        write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2881             iii=iii+1
2882             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2883               endif
2884           endif
2885         else if ((ions).and.(card(1:6).eq.'HETATM')) then
2886        if (firstion.eq.0) then
2887        firstion=1
2888        if (unres_pdb) then
2889          do j=1,3
2890            dc(j,ires)=sccor(j,iii)
2891          enddo
2892        else
2893           call sccenter(ires,iii,sccor)
2894        endif
2895        endif
2896           read (card(12:16),*) atom
2897           print *,"HETATOM", atom
2898           read (card(18:20),'(a3)') res
2899           if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
2900           (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
2901           .or.(atom(1:2).eq.'K ')) &
2902           then
2903            ires=ires+1
2904            if (molecule.ne.5) molecprev=molecule
2905            molecule=5
2906            nres_molec(molecule)=nres_molec(molecule)+1
2907            print *,"HERE",nres_molec(molecule)
2908            itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2909            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2910           endif
2911         endif !atom
2912       enddo
2913    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2914       if (ires.eq.0) return
2915 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2916 ! system
2917       nres=ires
2918       if ((ires_old.ne.ires).and.(molecule.ne.5)) &
2919          nres_molec(molecule)=nres_molec(molecule)-2
2920 !      print *,'I have', nres_molec(:)
2921
2922       do k=1,4 ! ions are without dummy
2923        if (nres_molec(k).eq.0) cycle
2924       do i=2,nres-1
2925 !        write (iout,*) i,itype(i,1)
2926 !        if (itype(i,1).eq.ntyp1) then
2927 !          write (iout,*) "dummy",i,itype(i,1)
2928 !          do j=1,3
2929 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2930 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
2931 !            dc(j,i)=c(j,i)
2932 !          enddo
2933 !        endif
2934         if (itype(i,k).eq.ntyp1_molec(k)) then
2935          if (itype(i+1,k).eq.ntyp1_molec(k)) then
2936           if (itype(i+2,k).eq.0) then
2937 !           print *,"masz sieczke"
2938            do j=1,5
2939             if (itype(i+2,j).ne.0) then
2940             itype(i+1,k)=0
2941             itype(i+1,j)=ntyp1_molec(j)
2942             nres_molec(k)=nres_molec(k)-1
2943             nres_molec(j)=nres_molec(j)+1
2944             go to 3331
2945             endif
2946            enddo
2947  3331      continue
2948           endif
2949 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2950 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
2951 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
2952            if (unres_pdb) then
2953 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2954 !            print *,i,'tu dochodze'
2955             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2956             if (fail) then
2957               e2(1)=0.0d0
2958               e2(2)=1.0d0
2959               e2(3)=0.0d0
2960             endif !fail
2961 !            print *,i,'a tu?'
2962             do j=1,3
2963              c(j,i)=c(j,i-1)-1.9d0*e2(j)
2964             enddo
2965            else   !unres_pdb
2966            do j=1,3
2967              dcj=(c(j,i-2)-c(j,i-3))/2.0
2968             if (dcj.eq.0) dcj=1.23591524223
2969              c(j,i)=c(j,i-1)+dcj
2970              c(j,nres+i)=c(j,i)
2971            enddo
2972           endif   !unres_pdb
2973          else     !itype(i+1,1).eq.ntyp1
2974           if (unres_pdb) then
2975 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2976             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2977             if (fail) then
2978               e2(1)=0.0d0
2979               e2(2)=1.0d0
2980               e2(3)=0.0d0
2981             endif
2982             do j=1,3
2983               c(j,i)=c(j,i+1)-1.9d0*e2(j)
2984             enddo
2985           else !unres_pdb
2986            do j=1,3
2987             dcj=(c(j,i+3)-c(j,i+2))/2.0
2988             if (dcj.eq.0) dcj=1.23591524223
2989             c(j,i)=c(j,i+1)-dcj
2990             c(j,nres+i)=c(j,i)
2991            enddo
2992           endif !unres_pdb
2993          endif !itype(i+1,1).eq.ntyp1
2994         endif  !itype.eq.ntyp1
2995
2996       enddo
2997       enddo
2998 ! Calculate the CM of the last side chain.
2999       if (iii.gt.0)  then
3000       if (unres_pdb) then
3001         do j=1,3
3002           dc(j,ires)=sccor(j,iii)
3003         enddo
3004       else
3005         call sccenter(ires,iii,sccor)
3006       endif
3007       endif
3008 !      nres=ires
3009       nsup=nres
3010       nstart_sup=1
3011 !      print *,"molecule",molecule
3012       if ((itype(nres,1).ne.10)) then
3013         nres=nres+1
3014           if (molecule.eq.5) molecule=molecprev
3015         itype(nres,molecule)=ntyp1_molec(molecule)
3016         nres_molec(molecule)=nres_molec(molecule)+1
3017         if (unres_pdb) then
3018 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3019           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3020           if (fail) then
3021             e2(1)=0.0d0
3022             e2(2)=1.0d0
3023             e2(3)=0.0d0
3024           endif
3025           do j=1,3
3026             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3027           enddo
3028         else
3029         do j=1,3
3030           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3031           c(j,nres)=c(j,nres-1)+dcj
3032           c(j,2*nres)=c(j,nres)
3033         enddo
3034         endif
3035       endif
3036 !     print *,'I have',nres, nres_molec(:)
3037
3038 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3039 #ifdef WHAM_RUN
3040       if (nres.ne.nres0) then
3041         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3042                        " NRES0=",nres0
3043         stop "Error nres value in WHAM input"
3044       endif
3045 #endif
3046 !---------------------------------
3047 !el reallocate tables
3048 !      do i=1,maxres/3
3049 !       do j=1,2
3050 !         hfrag_alloc(j,i)=hfrag(j,i)
3051 !        enddo
3052 !       do j=1,4
3053 !         bfrag_alloc(j,i)=bfrag(j,i)
3054 !        enddo
3055 !      enddo
3056
3057 !      deallocate(hfrag)
3058 !      deallocate(bfrag)
3059 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
3060 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3061 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3062 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
3063
3064 !      do i=1,nhfrag
3065 !       do j=1,2
3066 !         hfrag(j,i)=hfrag_alloc(j,i)
3067 !        enddo
3068 !      enddo
3069 !      do i=1,nbfrag
3070 !       do j=1,4
3071 !         bfrag(j,i)=bfrag_alloc(j,i)
3072 !        enddo
3073 !      enddo
3074 !el end reallocate tables
3075 !---------------------------------
3076       do i=2,nres-1
3077         do j=1,3
3078           c(j,i+nres)=dc(j,i)
3079         enddo
3080       enddo
3081       do j=1,3
3082         c(j,nres+1)=c(j,1)
3083         c(j,2*nres)=c(j,nres)
3084       enddo
3085
3086       if (itype(1,1).eq.ntyp1) then
3087         nsup=nsup-1
3088         nstart_sup=2
3089         if (unres_pdb) then
3090 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3091           call refsys(2,3,4,e1,e2,e3,fail)
3092           if (fail) then
3093             e2(1)=0.0d0
3094             e2(2)=1.0d0
3095             e2(3)=0.0d0
3096           endif
3097           do j=1,3
3098             c(j,1)=c(j,2)-1.9d0*e2(j)
3099           enddo
3100         else
3101         do j=1,3
3102           dcj=(c(j,4)-c(j,3))/2.0
3103           c(j,1)=c(j,2)-dcj
3104           c(j,nres+1)=c(j,1)
3105         enddo
3106         endif
3107       endif
3108 ! First lets assign correct dummy to correct type of chain
3109 ! 1) First residue
3110       if (itype(1,1).eq.ntyp1) then
3111         if (itype(2,1).eq.0) then
3112          do j=2,5
3113            if (itype(2,j).ne.0) then
3114            itype(1,1)=0
3115            itype(1,j)=ntyp1_molec(j)
3116            nres_molec(1)=nres_molec(1)-1
3117            nres_molec(j)=nres_molec(j)+1
3118            go to 3231
3119            endif
3120          enddo
3121 3231    continue
3122         endif
3123        endif
3124        print *,'I have',nres, nres_molec(:)
3125
3126 ! Copy the coordinates to reference coordinates
3127 !      do i=1,2*nres
3128 !        do j=1,3
3129 !          cref(j,i)=c(j,i)
3130 !        enddo
3131 !      enddo
3132 ! Calculate internal coordinates.
3133       if (lprn) then
3134       write (iout,'(/a)') &
3135         "Cartesian coordinates of the reference structure"
3136       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3137        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3138       do ires=1,nres
3139         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3140           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3141           (c(j,ires+nres),j=1,3)
3142       enddo
3143       endif
3144 ! znamy już nres wiec mozna alokowac tablice
3145 ! Calculate internal coordinates.
3146       if(me.eq.king.or..not.out1file)then
3147        write (iout,'(a)') &
3148          "Backbone and SC coordinates as read from the PDB"
3149        do ires=1,nres
3150         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3151           ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3152           (c(j,nres+ires),j=1,3)
3153        enddo
3154       endif
3155 ! NOW LETS ROCK! SORTING
3156       allocate(c_temporary(3,2*nres))
3157       allocate(itype_temporary(nres,5))
3158       allocate(molnum(nres))
3159       allocate(istype_temp(nres))
3160        itype_temporary(:,:)=0
3161       seqalingbegin=1
3162       do k=1,5
3163         do i=1,nres
3164          if (itype(i,k).ne.0) then
3165           do j=1,3
3166           c_temporary(j,seqalingbegin)=c(j,i)
3167           c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3168
3169           enddo
3170           itype_temporary(seqalingbegin,k)=itype(i,k)
3171           print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3172           istype_temp(seqalingbegin)=istype(i)
3173           molnum(seqalingbegin)=k
3174           seqalingbegin=seqalingbegin+1
3175          endif
3176         enddo
3177        enddo
3178        do i=1,2*nres
3179         do j=1,3
3180         c(j,i)=c_temporary(j,i)
3181         enddo
3182        enddo
3183        do k=1,5
3184         do i=1,nres
3185          itype(i,k)=itype_temporary(i,k)
3186          istype(i)=istype_temp(i)
3187         enddo
3188        enddo
3189       if (itype(1,1).eq.ntyp1) then
3190         nsup=nsup-1
3191         nstart_sup=2
3192         if (unres_pdb) then
3193 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3194           call refsys(2,3,4,e1,e2,e3,fail)
3195           if (fail) then
3196             e2(1)=0.0d0
3197             e2(2)=1.0d0
3198             e2(3)=0.0d0
3199           endif
3200           do j=1,3
3201             c(j,1)=c(j,2)-1.9d0*e2(j)
3202           enddo
3203         else
3204         do j=1,3
3205           dcj=(c(j,4)-c(j,3))/2.0
3206           c(j,1)=c(j,2)-dcj
3207           c(j,nres+1)=c(j,1)
3208         enddo
3209         endif
3210       endif
3211
3212       if (lprn) then
3213       write (iout,'(/a)') &
3214         "Cartesian coordinates of the reference structure after sorting"
3215       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3216        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3217       do ires=1,nres
3218         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3219           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3220           (c(j,ires+nres),j=1,3)
3221       enddo
3222       endif
3223
3224 !       print *,seqalingbegin,nres
3225       if(.not.allocated(vbld)) then
3226        allocate(vbld(2*nres))
3227        do i=1,2*nres
3228          vbld(i)=0.d0
3229        enddo
3230       endif
3231       if(.not.allocated(vbld_inv)) then
3232        allocate(vbld_inv(2*nres))
3233        do i=1,2*nres
3234          vbld_inv(i)=0.d0
3235        enddo
3236       endif
3237 !!!el
3238       if(.not.allocated(theta)) then
3239         allocate(theta(nres+2))
3240         theta(:)=0.0d0
3241       endif
3242
3243       if(.not.allocated(phi)) allocate(phi(nres+2))
3244       if(.not.allocated(alph)) allocate(alph(nres+2))
3245       if(.not.allocated(omeg)) allocate(omeg(nres+2))
3246       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3247       if(.not.allocated(phiref)) allocate(phiref(nres+2))
3248       if(.not.allocated(costtab)) allocate(costtab(nres))
3249       if(.not.allocated(sinttab)) allocate(sinttab(nres))
3250       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3251       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3252       if(.not.allocated(xxref)) allocate(xxref(nres))
3253       if(.not.allocated(yyref)) allocate(yyref(nres))
3254       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3255       if(.not.allocated(dc_norm)) then
3256 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3257         allocate(dc_norm(3,0:2*nres+2))
3258         dc_norm(:,:)=0.d0
3259       endif
3260
3261       call int_from_cart(.true.,.false.)
3262       call sc_loc_geom(.false.)
3263       do i=1,nres
3264         thetaref(i)=theta(i)
3265         phiref(i)=phi(i)
3266       enddo
3267 !      do i=1,2*nres
3268 !        vbld_inv(i)=0.d0
3269 !        vbld(i)=0.d0
3270 !      enddo
3271
3272       do i=1,nres-1
3273         do j=1,3
3274           dc(j,i)=c(j,i+1)-c(j,i)
3275           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3276         enddo
3277       enddo
3278       do i=2,nres-1
3279         do j=1,3
3280           dc(j,i+nres)=c(j,i+nres)-c(j,i)
3281           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3282         enddo
3283 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3284 !        vbld_inv(i+nres)
3285       enddo
3286 !      call chainbuild
3287 ! Copy the coordinates to reference coordinates
3288 ! Splits to single chain if occurs
3289
3290 !      do i=1,2*nres
3291 !        do j=1,3
3292 !          cref(j,i,cou)=c(j,i)
3293 !        enddo
3294 !      enddo
3295 !
3296       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3297       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3298       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3299 !-----------------------------
3300       kkk=1
3301       lll=0
3302       cou=1
3303         write (iout,*) "symetr", symetr
3304       do i=1,nres
3305       lll=lll+1
3306 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3307       if (i.gt.1) then
3308       if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3309       chain_length=lll-1
3310       kkk=kkk+1
3311 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3312       lll=1
3313       endif
3314       endif
3315         do j=1,3
3316           cref(j,i,cou)=c(j,i)
3317           cref(j,i+nres,cou)=c(j,i+nres)
3318           if (i.le.nres) then
3319           chain_rep(j,lll,kkk)=c(j,i)
3320           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3321           endif
3322          enddo
3323       enddo
3324       write (iout,*) chain_length
3325       if (chain_length.eq.0) chain_length=nres
3326       do j=1,3
3327       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3328       chain_rep(j,chain_length+nres,symetr) &
3329       =chain_rep(j,chain_length+nres,1)
3330       enddo
3331 ! diagnostic
3332 !       write (iout,*) "spraw lancuchy",chain_length,symetr
3333 !       do i=1,4
3334 !         do kkk=1,chain_length
3335 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3336 !         enddo
3337 !        enddo
3338 ! enddiagnostic
3339 ! makes copy of chains
3340         write (iout,*) "symetr", symetr
3341       do j=1,3
3342       dc(j,0)=c(j,1)
3343       enddo
3344
3345       if (symetr.gt.1) then
3346        call permut(symetr)
3347        nperm=1
3348        do i=1,symetr
3349        nperm=nperm*i
3350        enddo
3351        do i=1,nperm
3352        write(iout,*) (tabperm(i,kkk),kkk=1,4)
3353        enddo
3354        do i=1,nperm
3355         cou=0
3356         do kkk=1,symetr
3357          icha=tabperm(i,kkk)
3358 !         write (iout,*) i,icha
3359          do lll=1,chain_length
3360           cou=cou+1
3361            if (cou.le.nres) then
3362            do j=1,3
3363             kupa=mod(lll,chain_length)
3364             iprzes=(kkk-1)*chain_length+lll
3365             if (kupa.eq.0) kupa=chain_length
3366 !            write (iout,*) "kupa", kupa
3367             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3368             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3369           enddo
3370           endif
3371          enddo
3372         enddo
3373        enddo
3374        endif
3375 !-koniec robienia kopii
3376 ! diag
3377       do kkk=1,nperm
3378       write (iout,*) "nowa struktura", nperm
3379       do i=1,nres
3380         write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3381       cref(2,i,kkk),&
3382       cref(3,i,kkk),cref(1,nres+i,kkk),&
3383       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3384       enddo
3385   100 format (//'              alpha-carbon coordinates       ',&
3386                 '     centroid coordinates'/ &
3387                 '       ', 6X,'X',11X,'Y',11X,'Z', &
3388                                 10X,'X',11X,'Y',11X,'Z')
3389   110 format (a,'(',i3,')',6f12.5)
3390
3391       enddo
3392 !c enddiag
3393       do j=1,nbfrag
3394         do i=1,4
3395          bfrag(i,j)=bfrag(i,j)-ishift
3396         enddo
3397       enddo
3398
3399       do j=1,nhfrag
3400         do i=1,2
3401          hfrag(i,j)=hfrag(i,j)-ishift
3402         enddo
3403       enddo
3404       ishift_pdb=ishift
3405
3406       return
3407       end subroutine readpdb
3408 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3409 !-----------------------------------------------------------------------------
3410 ! readrtns_CSA.F
3411 !-----------------------------------------------------------------------------
3412       subroutine read_control
3413 !
3414 ! Read contorl data
3415 !
3416 !      use geometry_data
3417       use comm_machsw
3418       use energy_data
3419       use control_data
3420       use compare_data
3421       use MCM_data
3422       use map_data
3423       use csa_data
3424       use MD_data
3425       use MPI_data
3426       use random, only: random_init
3427 !      implicit real*8 (a-h,o-z)
3428 !      include 'DIMENSIONS'
3429 #ifdef MP
3430       use prng, only:prng_restart
3431       include 'mpif.h'
3432       logical :: OKRandom!, prng_restart
3433       real(kind=8) :: r1
3434 #endif
3435 !      include 'COMMON.IOUNITS'
3436 !      include 'COMMON.TIME1'
3437 !      include 'COMMON.THREAD'
3438 !      include 'COMMON.SBRIDGE'
3439 !      include 'COMMON.CONTROL'
3440 !      include 'COMMON.MCM'
3441 !      include 'COMMON.MAP'
3442 !      include 'COMMON.HEADER'
3443 !      include 'COMMON.CSA'
3444 !      include 'COMMON.CHAIN'
3445 !      include 'COMMON.MUCA'
3446 !      include 'COMMON.MD'
3447 !      include 'COMMON.FFIELD'
3448 !      include 'COMMON.INTERACT'
3449 !      include 'COMMON.SETUP'
3450 !el      integer :: KDIAG,ICORFL,IXDR
3451 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3452       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3453         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
3454 !      character(len=80) :: ucase
3455       character(len=640) :: controlcard
3456
3457       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3458       integer i
3459
3460       nglob_csa=0
3461       eglob_csa=1d99
3462       nmin_csa=0
3463       read (INP,'(a)') titel
3464       call card_concat(controlcard,.true.)
3465 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3466 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3467       call reada(controlcard,'SEED',seed,0.0D0)
3468       call random_init(seed)
3469 ! Set up the time limit (caution! The time must be input in minutes!)
3470       read_cart=index(controlcard,'READ_CART').gt.0
3471       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3472       call readi(controlcard,'SYM',symetr,1)
3473       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3474       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3475       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3476       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3477       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3478       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3479       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3480       call reada(controlcard,'DRMS',drms,0.1D0)
3481       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3482        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3483        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3484        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3485        write (iout,'(a,f10.1)')'DRMS    = ',drms
3486        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3487        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3488       endif
3489       call readi(controlcard,'NZ_START',nz_start,0)
3490       call readi(controlcard,'NZ_END',nz_end,0)
3491       call readi(controlcard,'IZ_SC',iz_sc,0)
3492       timlim=60.0D0*timlim
3493       safety = 60.0d0*safety
3494       timem=timlim
3495       modecalc=0
3496       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3497 !C SHIELD keyword sets if the shielding effect of side-chains is used
3498 !C 0 denots no shielding is used all peptide are equally despite the
3499 !C solvent accesible area
3500 !C 1 the newly introduced function
3501 !C 2 reseved for further possible developement
3502       call readi(controlcard,'SHIELD',shield_mode,0)
3503 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3504         write(iout,*) "shield_mode",shield_mode
3505 !C  Varibles set size of box
3506       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3507       protein=index(controlcard,"PROTEIN").gt.0
3508       ions=index(controlcard,"IONS").gt.0
3509       nucleic=index(controlcard,"NUCLEIC").gt.0
3510       write (iout,*) "with_theta_constr ",with_theta_constr
3511       AFMlog=(index(controlcard,'AFM'))
3512       selfguide=(index(controlcard,'SELFGUIDE'))
3513       print *,'AFMlog',AFMlog,selfguide,"KUPA"
3514       call readi(controlcard,'GENCONSTR',genconstr,0)
3515       call reada(controlcard,'BOXX',boxxsize,100.0d0)
3516       call reada(controlcard,'BOXY',boxysize,100.0d0)
3517       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3518       call readi(controlcard,'TUBEMOD',tubemode,0)
3519       write (iout,*) TUBEmode,"TUBEMODE"
3520       if (TUBEmode.gt.0) then
3521        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3522        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3523        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3524        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3525        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3526        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3527        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3528        buftubebot=bordtubebot+tubebufthick
3529        buftubetop=bordtubetop-tubebufthick
3530       endif
3531
3532 ! CUTOFFF ON ELECTROSTATICS
3533       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3534       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3535       write(iout,*) "R_CUT_ELE=",r_cut_ele
3536 ! Lipidic parameters
3537       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3538       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3539       if (lipthick.gt.0.0d0) then
3540        bordliptop=(boxzsize+lipthick)/2.0
3541        bordlipbot=bordliptop-lipthick
3542       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3543       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3544       buflipbot=bordlipbot+lipbufthick
3545       bufliptop=bordliptop-lipbufthick
3546       if ((lipbufthick*2.0d0).gt.lipthick) &
3547        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3548       endif !lipthick.gt.0
3549       write(iout,*) "bordliptop=",bordliptop
3550       write(iout,*) "bordlipbot=",bordlipbot
3551       write(iout,*) "bufliptop=",bufliptop
3552       write(iout,*) "buflipbot=",buflipbot
3553       write (iout,*) "SHIELD MODE",shield_mode
3554
3555 !C-------------------------
3556       minim=(index(controlcard,'MINIMIZE').gt.0)
3557       dccart=(index(controlcard,'CART').gt.0)
3558       overlapsc=(index(controlcard,'OVERLAP').gt.0)
3559       overlapsc=.not.overlapsc
3560       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3561       searchsc=.not.searchsc
3562       sideadd=(index(controlcard,'SIDEADD').gt.0)
3563       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3564       outpdb=(index(controlcard,'PDBOUT').gt.0)
3565       outmol2=(index(controlcard,'MOL2OUT').gt.0)
3566       pdbref=(index(controlcard,'PDBREF').gt.0)
3567       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3568       indpdb=index(controlcard,'PDBSTART')
3569       extconf=(index(controlcard,'EXTCONF').gt.0)
3570       call readi(controlcard,'IPRINT',iprint,0)
3571       call readi(controlcard,'MAXGEN',maxgen,10000)
3572       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3573       call readi(controlcard,"KDIAG",kdiag,0)
3574       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3575       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3576        write (iout,*) "RESCALE_MODE",rescale_mode
3577       split_ene=index(controlcard,'SPLIT_ENE').gt.0
3578       if (index(controlcard,'REGULAR').gt.0.0D0) then
3579         call reada(controlcard,'WEIDIS',weidis,0.1D0)
3580         modecalc=1
3581         refstr=.true.
3582       endif
3583       if (index(controlcard,'CHECKGRAD').gt.0) then
3584         modecalc=5
3585         if (index(controlcard,'CART').gt.0) then
3586           icheckgrad=1
3587         elseif (index(controlcard,'CARINT').gt.0) then
3588           icheckgrad=2
3589         else
3590           icheckgrad=3
3591         endif
3592       elseif (index(controlcard,'THREAD').gt.0) then
3593         modecalc=2
3594         call readi(controlcard,'THREAD',nthread,0)
3595         if (nthread.gt.0) then
3596           call reada(controlcard,'WEIDIS',weidis,0.1D0)
3597         else
3598           if (fg_rank.eq.0) &
3599           write (iout,'(a)')'A number has to follow the THREAD keyword.'
3600           stop 'Error termination in Read_Control.'
3601         endif
3602       else if (index(controlcard,'MCMA').gt.0) then
3603         modecalc=3
3604       else if (index(controlcard,'MCEE').gt.0) then
3605         modecalc=6
3606       else if (index(controlcard,'MULTCONF').gt.0) then
3607         modecalc=4
3608       else if (index(controlcard,'MAP').gt.0) then
3609         modecalc=7
3610         call readi(controlcard,'MAP',nmap,0)
3611       else if (index(controlcard,'CSA').gt.0) then
3612         modecalc=8
3613 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3614 !rc
3615 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3616 !rc
3617 !rc        modecalc=9
3618 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3619 !fmc        modecalc=10
3620       else if (index(controlcard,'SOFTREG').gt.0) then
3621         modecalc=11
3622       else if (index(controlcard,'CHECK_BOND').gt.0) then
3623         modecalc=-1
3624       else if (index(controlcard,'TEST').gt.0) then
3625         modecalc=-2
3626       else if (index(controlcard,'MD').gt.0) then
3627         modecalc=12
3628       else if (index(controlcard,'RE ').gt.0) then
3629         modecalc=14
3630       endif
3631
3632       lmuca=index(controlcard,'MUCA').gt.0
3633       call readi(controlcard,'MUCADYN',mucadyn,0)
3634       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3635       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3636        then
3637        write (iout,*) 'MUCADYN=',mucadyn
3638        write (iout,*) 'MUCASMOOTH=',muca_smooth
3639       endif
3640
3641       iscode=index(controlcard,'ONE_LETTER')
3642       indphi=index(controlcard,'PHI')
3643       indback=index(controlcard,'BACK')
3644       iranconf=index(controlcard,'RAND_CONF')
3645       i2ndstr=index(controlcard,'USE_SEC_PRED')
3646       gradout=index(controlcard,'GRADOUT').gt.0
3647       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3648       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3649       if (me.eq.king .or. .not.out1file ) &
3650         write (iout,*) "DISTCHAINMAX",distchainmax
3651
3652       if(me.eq.king.or..not.out1file) &
3653        write (iout,'(2a)') diagmeth(kdiag),&
3654         ' routine used to diagonalize matrices.'
3655       if (shield_mode.gt.0) then
3656       pi=3.141592d0
3657 !C VSolvSphere the volume of solving sphere
3658 !C      print *,pi,"pi"
3659 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3660 !C there will be no distinction between proline peptide group and normal peptide
3661 !C group in case of shielding parameters
3662       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3663       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3664       write (iout,*) VSolvSphere,VSolvSphere_div
3665 !C long axis of side chain
3666       do i=1,ntyp
3667       long_r_sidechain(i)=vbldsc0(1,i)
3668       short_r_sidechain(i)=sigma0(i)
3669       write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3670          sigma0(i)
3671       enddo
3672       buff_shield=1.0d0
3673       endif
3674       return
3675       end subroutine read_control
3676 !-----------------------------------------------------------------------------
3677       subroutine read_REMDpar
3678 !
3679 ! Read REMD settings
3680 !
3681 !       use control
3682 !       use energy
3683 !       use geometry
3684       use REMD_data
3685       use MPI_data
3686       use control_data, only:out1file
3687 !      implicit real*8 (a-h,o-z)
3688 !      include 'DIMENSIONS'
3689 !      include 'COMMON.IOUNITS'
3690 !      include 'COMMON.TIME1'
3691 !      include 'COMMON.MD'
3692       use MD_data
3693 !el #ifndef LANG0
3694 !el      include 'COMMON.LANGEVIN'
3695 !el #else
3696 !el      include 'COMMON.LANGEVIN.lang0'
3697 !el #endif
3698 !      include 'COMMON.INTERACT'
3699 !      include 'COMMON.NAMES'
3700 !      include 'COMMON.GEO'
3701 !      include 'COMMON.REMD'
3702 !      include 'COMMON.CONTROL'
3703 !      include 'COMMON.SETUP'
3704 !      character(len=80) :: ucase
3705       character(len=320) :: controlcard
3706       character(len=3200) :: controlcard1
3707       integer :: iremd_m_total
3708 !el local variables
3709       integer :: i
3710 !     real(kind=8) :: var,ene
3711
3712       if(me.eq.king.or..not.out1file) &
3713        write (iout,*) "REMD setup"
3714
3715       call card_concat(controlcard,.true.)
3716       call readi(controlcard,"NREP",nrep,3)
3717       call readi(controlcard,"NSTEX",nstex,1000)
3718       call reada(controlcard,"RETMIN",retmin,10.0d0)
3719       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3720       mremdsync=(index(controlcard,'SYNC').gt.0)
3721       call readi(controlcard,"NSYN",i_sync_step,100)
3722       restart1file=(index(controlcard,'REST1FILE').gt.0)
3723       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3724       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3725       if(max_cache_traj_use.gt.max_cache_traj) &
3726                 max_cache_traj_use=max_cache_traj
3727       if(me.eq.king.or..not.out1file) then
3728 !d       if (traj1file) then
3729 !rc caching is in testing - NTWX is not ignored
3730 !d        write (iout,*) "NTWX value is ignored"
3731 !d        write (iout,*) "  trajectory is stored to one file by master"
3732 !d        write (iout,*) "  before exchange at NSTEX intervals"
3733 !d       endif
3734        write (iout,*) "NREP= ",nrep
3735        write (iout,*) "NSTEX= ",nstex
3736        write (iout,*) "SYNC= ",mremdsync
3737        write (iout,*) "NSYN= ",i_sync_step
3738        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3739       endif
3740       remd_tlist=.false.
3741       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3742       if (index(controlcard,'TLIST').gt.0) then
3743          remd_tlist=.true.
3744          call card_concat(controlcard1,.true.)
3745          read(controlcard1,*) (remd_t(i),i=1,nrep)
3746          if(me.eq.king.or..not.out1file) &
3747           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3748       endif
3749       remd_mlist=.false.
3750       if (index(controlcard,'MLIST').gt.0) then
3751          remd_mlist=.true.
3752          call card_concat(controlcard1,.true.)
3753          read(controlcard1,*) (remd_m(i),i=1,nrep)
3754          if(me.eq.king.or..not.out1file) then
3755           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3756           iremd_m_total=0
3757           do i=1,nrep
3758            iremd_m_total=iremd_m_total+remd_m(i)
3759           enddo
3760           write (iout,*) 'Total number of replicas ',iremd_m_total
3761          endif
3762       endif
3763       if(me.eq.king.or..not.out1file) &
3764        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3765       return
3766       end subroutine read_REMDpar
3767 !-----------------------------------------------------------------------------
3768       subroutine read_MDpar
3769 !
3770 ! Read MD settings
3771 !
3772       use control_data, only: r_cut,rlamb,out1file
3773       use energy_data
3774       use geometry_data, only: pi
3775       use MPI_data
3776 !      implicit real*8 (a-h,o-z)
3777 !      include 'DIMENSIONS'
3778 !      include 'COMMON.IOUNITS'
3779 !      include 'COMMON.TIME1'
3780 !      include 'COMMON.MD'
3781       use MD_data
3782 !el #ifndef LANG0
3783 !el      include 'COMMON.LANGEVIN'
3784 !el #else
3785 !el      include 'COMMON.LANGEVIN.lang0'
3786 !el #endif
3787 !      include 'COMMON.INTERACT'
3788 !      include 'COMMON.NAMES'
3789 !      include 'COMMON.GEO'
3790 !      include 'COMMON.SETUP'
3791 !      include 'COMMON.CONTROL'
3792 !      include 'COMMON.SPLITELE'
3793 !      character(len=80) :: ucase
3794       character(len=320) :: controlcard
3795 !el local variables
3796       integer :: i,j
3797       real(kind=8) :: eta
3798
3799       call card_concat(controlcard,.true.)
3800       call readi(controlcard,"NSTEP",n_timestep,1000000)
3801       call readi(controlcard,"NTWE",ntwe,100)
3802       call readi(controlcard,"NTWX",ntwx,1000)
3803       call reada(controlcard,"DT",d_time,1.0d-1)
3804       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3805       call reada(controlcard,"DAMAX",damax,1.0d1)
3806       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3807       call readi(controlcard,"LANG",lang,0)
3808       RESPA = index(controlcard,"RESPA") .gt. 0
3809       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3810       ntime_split0=ntime_split
3811       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3812       ntime_split0=ntime_split
3813       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3814       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3815       rest = index(controlcard,"REST").gt.0
3816       tbf = index(controlcard,"TBF").gt.0
3817       usampl = index(controlcard,"USAMPL").gt.0
3818       mdpdb = index(controlcard,"MDPDB").gt.0
3819       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3820       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3821       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3822       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3823       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3824       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3825       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3826       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3827       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3828       large = index(controlcard,"LARGE").gt.0
3829       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3830       rattle = index(controlcard,"RATTLE").gt.0
3831 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3832       nset=0
3833       if(usampl) then
3834         call read_fragments
3835       endif
3836
3837       if(me.eq.king.or..not.out1file) then
3838        write (iout,*)
3839        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3840        write (iout,*)
3841        write (iout,'(a)') "The units are:"
3842        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3843        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3844         " acceleration: angstrom/(48.9 fs)**2"
3845        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3846        write (iout,*)
3847        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3848        write (iout,'(a60,f10.5,a)') &
3849         "Initial time step of numerical integration:",d_time,&
3850         " natural units"
3851        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3852        if (RESPA) then
3853         write (iout,'(2a,i4,a)') &
3854           "A-MTS algorithm used; initial time step for fast-varying",&
3855           " short-range forces split into",ntime_split," steps."
3856         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3857          r_cut," lambda",rlamb
3858        endif
3859        write (iout,'(2a,f10.5)') &
3860         "Maximum acceleration threshold to reduce the time step",&
3861         "/increase split number:",damax
3862        write (iout,'(2a,f10.5)') &
3863         "Maximum predicted energy drift to reduce the timestep",&
3864         "/increase split number:",edriftmax
3865        write (iout,'(a60,f10.5)') &
3866        "Maximum velocity threshold to reduce velocities:",dvmax
3867        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3868        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3869        if (rattle) write (iout,'(a60)') &
3870         "Rattle algorithm used to constrain the virtual bonds"
3871       endif
3872       reset_fricmat=1000
3873       if (lang.gt.0) then
3874         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3875         call reada(controlcard,"RWAT",rwat,1.4d0)
3876         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3877         surfarea=index(controlcard,"SURFAREA").gt.0
3878         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3879         if(me.eq.king.or..not.out1file)then
3880          write (iout,'(/a,$)') "Langevin dynamics calculation"
3881          if (lang.eq.1) then
3882           write (iout,'(a/)') &
3883             " with direct integration of Langevin equations"
3884          else if (lang.eq.2) then
3885           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3886          else if (lang.eq.3) then
3887           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3888          else if (lang.eq.4) then
3889           write (iout,'(a/)') " in overdamped mode"
3890          else
3891           write (iout,'(//a,i5)') &
3892             "=========== ERROR: Unknown Langevin dynamics mode:",lang
3893           stop
3894          endif
3895          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3896          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3897          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3898          write (iout,'(a60,f10.5)') &
3899          "Scaling factor of the friction forces:",scal_fric
3900          if (surfarea) write (iout,'(2a,i10,a)') &
3901            "Friction coefficients will be scaled by solvent-accessible",&
3902            " surface area every",reset_fricmat," steps."
3903         endif
3904 ! Calculate friction coefficients and bounds of stochastic forces
3905         eta=6*pi*cPoise*etawat
3906         if(me.eq.king.or..not.out1file) &
3907          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3908           eta
3909 !        allocate(gamp
3910         do j=1,5 !types of molecules
3911         gamp(j)=scal_fric*(pstok(j)+rwat)*eta
3912         stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
3913         enddo
3914         allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
3915         do j=1,5 !types of molecules
3916         do i=1,ntyp
3917           gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
3918           stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
3919         enddo
3920         enddo
3921
3922         if(me.eq.king.or..not.out1file)then
3923          write (iout,'(/2a/)') &
3924          "Radii of site types and friction coefficients and std's of",&
3925          " stochastic forces of fully exposed sites"
3926          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
3927          do i=1,ntyp
3928           write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
3929            gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
3930          enddo
3931         endif
3932       else if (tbf) then
3933         if(me.eq.king.or..not.out1file)then
3934          write (iout,'(a)') "Berendsen bath calculation"
3935          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3936          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3937          if (reset_moment) &
3938          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3939          count_reset_moment," steps"
3940          if (reset_vel) &
3941           write (iout,'(a,i10,a)') &
3942           "Velocities will be reset at random every",count_reset_vel,&
3943          " steps"
3944         endif
3945       else
3946         if(me.eq.king.or..not.out1file) &
3947          write (iout,'(a31)') "Microcanonical mode calculation"
3948       endif
3949       if(me.eq.king.or..not.out1file)then
3950        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3951        if (usampl) then
3952           write(iout,*) "MD running with constraints."
3953           write(iout,*) "Equilibration time ", eq_time, " mtus."
3954           write(iout,*) "Constraining ", nfrag," fragments."
3955           write(iout,*) "Length of each fragment, weight and q0:"
3956           do iset=1,nset
3957            write (iout,*) "Set of restraints #",iset
3958            do i=1,nfrag
3959               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3960                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3961            enddo
3962            write(iout,*) "constraints between ", npair, "fragments."
3963            write(iout,*) "constraint pairs, weights and q0:"
3964            do i=1,npair
3965             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3966                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3967            enddo
3968            write(iout,*) "angle constraints within ", nfrag_back,&
3969             "backbone fragments."
3970            write(iout,*) "fragment, weights:"
3971            do i=1,nfrag_back
3972             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3973                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3974                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3975            enddo
3976           enddo
3977         iset=mod(kolor,nset)+1
3978        endif
3979       endif
3980       if(me.eq.king.or..not.out1file) &
3981        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3982       return
3983       end subroutine read_MDpar
3984 !-----------------------------------------------------------------------------
3985       subroutine map_read
3986
3987       use map_data
3988 !      implicit real*8 (a-h,o-z)
3989 !      include 'DIMENSIONS'
3990 !      include 'COMMON.MAP'
3991 !      include 'COMMON.IOUNITS'
3992       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3993       character(len=80) :: mapcard      !,ucase
3994 !el local variables
3995       integer :: imap
3996 !     real(kind=8) :: var,ene
3997
3998       do imap=1,nmap
3999         read (inp,'(a)') mapcard
4000         mapcard=ucase(mapcard)
4001         if (index(mapcard,'PHI').gt.0) then
4002           kang(imap)=1
4003         else if (index(mapcard,'THE').gt.0) then
4004           kang(imap)=2
4005         else if (index(mapcard,'ALP').gt.0) then
4006           kang(imap)=3
4007         else if (index(mapcard,'OME').gt.0) then
4008           kang(imap)=4
4009         else
4010           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4011           stop 'Error - illegal variable spec in MAP card.'
4012         endif
4013         call readi (mapcard,'RES1',res1(imap),0)
4014         call readi (mapcard,'RES2',res2(imap),0)
4015         if (res1(imap).eq.0) then
4016           res1(imap)=res2(imap)
4017         else if (res2(imap).eq.0) then
4018           res2(imap)=res1(imap)
4019         endif
4020         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4021           write (iout,'(a)') &
4022           'Error - illegal definition of variable group in MAP.'
4023           stop 'Error - illegal definition of variable group in MAP.'
4024         endif
4025         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4026         call reada(mapcard,'TO',ang_to(imap),0.0D0)
4027         call readi(mapcard,'NSTEP',nstep(imap),0)
4028         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4029           write (iout,'(a)') &
4030            'Illegal boundary and/or step size specification in MAP.'
4031           stop 'Illegal boundary and/or step size specification in MAP.'
4032         endif
4033       enddo ! imap
4034       return
4035       end subroutine map_read
4036 !-----------------------------------------------------------------------------
4037       subroutine csaread
4038
4039       use control_data, only: vdisulf
4040       use csa_data
4041 !      implicit real*8 (a-h,o-z)
4042 !      include 'DIMENSIONS'
4043 !      include 'COMMON.IOUNITS'
4044 !      include 'COMMON.GEO'
4045 !      include 'COMMON.CSA'
4046 !      include 'COMMON.BANK'
4047 !      include 'COMMON.CONTROL'
4048 !      character(len=80) :: ucase
4049       character(len=620) :: mcmcard
4050 !el local variables
4051 !     integer :: ntf,ik,iw_pdb
4052 !     real(kind=8) :: var,ene
4053
4054       call card_concat(mcmcard,.true.)
4055
4056       call readi(mcmcard,'NCONF',nconf,50)
4057       call readi(mcmcard,'NADD',nadd,0)
4058       call readi(mcmcard,'JSTART',jstart,1)
4059       call readi(mcmcard,'JEND',jend,1)
4060       call readi(mcmcard,'NSTMAX',nstmax,500000)
4061       call readi(mcmcard,'N0',n0,1)
4062       call readi(mcmcard,'N1',n1,6)
4063       call readi(mcmcard,'N2',n2,4)
4064       call readi(mcmcard,'N3',n3,0)
4065       call readi(mcmcard,'N4',n4,0)
4066       call readi(mcmcard,'N5',n5,0)
4067       call readi(mcmcard,'N6',n6,10)
4068       call readi(mcmcard,'N7',n7,0)
4069       call readi(mcmcard,'N8',n8,0)
4070       call readi(mcmcard,'N9',n9,0)
4071       call readi(mcmcard,'N14',n14,0)
4072       call readi(mcmcard,'N15',n15,0)
4073       call readi(mcmcard,'N16',n16,0)
4074       call readi(mcmcard,'N17',n17,0)
4075       call readi(mcmcard,'N18',n18,0)
4076
4077       vdisulf=(index(mcmcard,'DYNSS').gt.0)
4078
4079       call readi(mcmcard,'NDIFF',ndiff,2)
4080       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4081       call readi(mcmcard,'IS1',is1,1)
4082       call readi(mcmcard,'IS2',is2,8)
4083       call readi(mcmcard,'NRAN0',nran0,4)
4084       call readi(mcmcard,'NRAN1',nran1,2)
4085       call readi(mcmcard,'IRR',irr,1)
4086       call readi(mcmcard,'NSEED',nseed,20)
4087       call readi(mcmcard,'NTOTAL',ntotal,10000)
4088       call reada(mcmcard,'CUT1',cut1,2.0d0)
4089       call reada(mcmcard,'CUT2',cut2,5.0d0)
4090       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4091       call readi(mcmcard,'ICMAX',icmax,3)
4092       call readi(mcmcard,'IRESTART',irestart,0)
4093 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
4094       ntbankm=0
4095 !!bankt
4096       call reada(mcmcard,'DELE',dele,20.0d0)
4097       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4098       call readi(mcmcard,'IREF',iref,0)
4099       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4100       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4101       call readi(mcmcard,'NCONF_IN',nconf_in,0)
4102       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4103       write (iout,*) "NCONF_IN",nconf_in
4104       return
4105       end subroutine csaread
4106 !-----------------------------------------------------------------------------
4107       subroutine mcmread
4108
4109       use mcm_data
4110       use control_data, only: MaxMoveType
4111       use MD_data
4112       use minim_data
4113 !      implicit real*8 (a-h,o-z)
4114 !      include 'DIMENSIONS'
4115 !      include 'COMMON.MCM'
4116 !      include 'COMMON.MCE'
4117 !      include 'COMMON.IOUNITS'
4118 !      character(len=80) :: ucase
4119       character(len=320) :: mcmcard
4120 !el local variables
4121       integer :: i
4122 !     real(kind=8) :: var,ene
4123
4124       call card_concat(mcmcard,.true.)
4125       call readi(mcmcard,'MAXACC',maxacc,100)
4126       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4127       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4128       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4129       call readi(mcmcard,'MAXREPM',maxrepm,200)
4130       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4131       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4132       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4133       call reada(mcmcard,'E_UP',e_up,5.0D0)
4134       call reada(mcmcard,'DELTE',delte,0.1D0)
4135       call readi(mcmcard,'NSWEEP',nsweep,5)
4136       call readi(mcmcard,'NSTEPH',nsteph,0)
4137       call readi(mcmcard,'NSTEPC',nstepc,0)
4138       call reada(mcmcard,'TMIN',tmin,298.0D0)
4139       call reada(mcmcard,'TMAX',tmax,298.0D0)
4140       call readi(mcmcard,'NWINDOW',nwindow,0)
4141       call readi(mcmcard,'PRINT_MC',print_mc,0)
4142       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4143       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4144       ent_read=(index(mcmcard,'ENT_READ').gt.0)
4145       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4146       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4147       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4148       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4149       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4150       if (nwindow.gt.0) then
4151         allocate(winstart(nwindow))     !!el (maxres)
4152         allocate(winend(nwindow))       !!el
4153         allocate(winlen(nwindow))       !!el
4154         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4155         do i=1,nwindow
4156           winlen(i)=winend(i)-winstart(i)+1
4157         enddo
4158       endif
4159       if (tmax.lt.tmin) tmax=tmin
4160       if (tmax.eq.tmin) then
4161         nstepc=0
4162         nsteph=0
4163       endif
4164       if (nstepc.gt.0 .and. nsteph.gt.0) then
4165         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4166         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4167       endif
4168       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4169 ! Probabilities of different move types
4170       sumpro_type(0)=0.0D0
4171       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4172       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4173       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4174       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
4175       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4176       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4177       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4178       do i=1,MaxMoveType
4179         print *,'i',i,' sumprotype',sumpro_type(i)
4180         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4181         print *,'i',i,' sumprotype',sumpro_type(i)
4182       enddo
4183       return
4184       end subroutine mcmread
4185 !-----------------------------------------------------------------------------
4186       subroutine read_minim
4187
4188       use minim_data
4189 !      implicit real*8 (a-h,o-z)
4190 !      include 'DIMENSIONS'
4191 !      include 'COMMON.MINIM'
4192 !      include 'COMMON.IOUNITS'
4193 !      character(len=80) :: ucase
4194       character(len=320) :: minimcard
4195 !el local variables
4196 !     integer :: ntf,ik,iw_pdb
4197 !     real(kind=8) :: var,ene
4198
4199       call card_concat(minimcard,.true.)
4200       call readi(minimcard,'MAXMIN',maxmin,2000)
4201       call readi(minimcard,'MAXFUN',maxfun,5000)
4202       call readi(minimcard,'MINMIN',minmin,maxmin)
4203       call readi(minimcard,'MINFUN',minfun,maxmin)
4204       call reada(minimcard,'TOLF',tolf,1.0D-2)
4205       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4206       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4207       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4208       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4209       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4210                'Options in energy minimization:'
4211       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4212        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4213        'MinMin:',MinMin,' MinFun:',MinFun,&
4214        ' TolF:',TolF,' RTolF:',RTolF
4215       return
4216       end subroutine read_minim
4217 !-----------------------------------------------------------------------------
4218       subroutine openunits
4219
4220       use MD_data, only: usampl
4221       use csa_data
4222       use MPI_data
4223       use control_data, only:out1file
4224       use control, only: getenv_loc
4225 !      implicit real*8 (a-h,o-z)
4226 !      include 'DIMENSIONS'
4227 #ifdef MPI
4228       include 'mpif.h'
4229       character(len=16) :: form,nodename
4230       integer :: nodelen,ierror,npos
4231 #endif
4232 !      include 'COMMON.SETUP'
4233 !      include 'COMMON.IOUNITS'
4234 !      include 'COMMON.MD'
4235 !      include 'COMMON.CONTROL'
4236       integer :: lenpre,lenpot,lentmp   !,ilen
4237 !el      external ilen
4238       character(len=3) :: out1file_text !,ucase
4239       character(len=3) :: ll
4240 !el      external ucase
4241 !el local variables
4242 !     integer :: ntf,ik,iw_pdb
4243 !     real(kind=8) :: var,ene
4244 !
4245 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4246       call getenv_loc("PREFIX",prefix)
4247       pref_orig = prefix
4248       call getenv_loc("POT",pot)
4249       call getenv_loc("DIRTMP",tmpdir)
4250       call getenv_loc("CURDIR",curdir)
4251       call getenv_loc("OUT1FILE",out1file_text)
4252 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4253       out1file_text=ucase(out1file_text)
4254       if (out1file_text(1:1).eq."Y") then
4255         out1file=.true.
4256       else
4257         out1file=fg_rank.gt.0
4258       endif
4259       lenpre=ilen(prefix)
4260       lenpot=ilen(pot)
4261       lentmp=ilen(tmpdir)
4262       if (lentmp.gt.0) then
4263           write (*,'(80(1h!))')
4264           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
4265           write (*,'(80(1h!))')
4266           write (*,*)"All output files will be on node /tmp directory."
4267 #ifdef MPI
4268         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4269         if (me.eq.king) then
4270           write (*,*) "The master node is ",nodename
4271         else if (fg_rank.eq.0) then
4272           write (*,*) "I am the CG slave node ",nodename
4273         else
4274           write (*,*) "I am the FG slave node ",nodename
4275         endif
4276 #endif
4277         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4278         lenpre = lentmp+lenpre+1
4279       endif
4280       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4281 ! Get the names and open the input files
4282 #if defined(WINIFL) || defined(WINPGI)
4283       open(1,file=pref_orig(:ilen(pref_orig))// &
4284         '.inp',status='old',readonly,shared)
4285        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4286 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4287 ! Get parameter filenames and open the parameter files.
4288       call getenv_loc('BONDPAR',bondname)
4289       open (ibond,file=bondname,status='old',readonly,shared)
4290       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4291       open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4292       call getenv_loc('THETPAR',thetname)
4293       open (ithep,file=thetname,status='old',readonly,shared)
4294       call getenv_loc('ROTPAR',rotname)
4295       open (irotam,file=rotname,status='old',readonly,shared)
4296       call getenv_loc('TORPAR',torname)
4297       open (itorp,file=torname,status='old',readonly,shared)
4298       call getenv_loc('TORDPAR',tordname)
4299       open (itordp,file=tordname,status='old',readonly,shared)
4300       call getenv_loc('FOURIER',fouriername)
4301       open (ifourier,file=fouriername,status='old',readonly,shared)
4302       call getenv_loc('ELEPAR',elename)
4303       open (ielep,file=elename,status='old',readonly,shared)
4304       call getenv_loc('SIDEPAR',sidename)
4305       open (isidep,file=sidename,status='old',readonly,shared)
4306
4307       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4308       open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4309       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4310       open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4311       call getenv_loc('TORPAR_NUCL',torname_nucl)
4312       open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4313       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4314       open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4315       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4316       open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4317
4318
4319 #elif (defined CRAY) || (defined AIX)
4320       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4321         action='read')
4322 !      print *,"Processor",myrank," opened file 1"
4323       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4324 !      print *,"Processor",myrank," opened file 9"
4325 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4326 ! Get parameter filenames and open the parameter files.
4327       call getenv_loc('BONDPAR',bondname)
4328       open (ibond,file=bondname,status='old',action='read')
4329       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4330       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4331
4332 !      print *,"Processor",myrank," opened file IBOND"
4333       call getenv_loc('THETPAR',thetname)
4334       open (ithep,file=thetname,status='old',action='read')
4335 !      print *,"Processor",myrank," opened file ITHEP"
4336       call getenv_loc('ROTPAR',rotname)
4337       open (irotam,file=rotname,status='old',action='read')
4338 !      print *,"Processor",myrank," opened file IROTAM"
4339       call getenv_loc('TORPAR',torname)
4340       open (itorp,file=torname,status='old',action='read')
4341 !      print *,"Processor",myrank," opened file ITORP"
4342       call getenv_loc('TORDPAR',tordname)
4343       open (itordp,file=tordname,status='old',action='read')
4344 !      print *,"Processor",myrank," opened file ITORDP"
4345       call getenv_loc('SCCORPAR',sccorname)
4346       open (isccor,file=sccorname,status='old',action='read')
4347 !      print *,"Processor",myrank," opened file ISCCOR"
4348       call getenv_loc('FOURIER',fouriername)
4349       open (ifourier,file=fouriername,status='old',action='read')
4350 !      print *,"Processor",myrank," opened file IFOURIER"
4351       call getenv_loc('ELEPAR',elename)
4352       open (ielep,file=elename,status='old',action='read')
4353 !      print *,"Processor",myrank," opened file IELEP"
4354       call getenv_loc('SIDEPAR',sidename)
4355       open (isidep,file=sidename,status='old',action='read')
4356
4357       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4358       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4359       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4360       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4361       call getenv_loc('TORPAR_NUCL',torname_nucl)
4362       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4363       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4364       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4365       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4366       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4367
4368       call getenv_loc('LIPTRANPAR',liptranname)
4369       open (iliptranpar,file=liptranname,status='old',action='read')
4370       call getenv_loc('TUBEPAR',tubename)
4371       open (itube,file=tubename,status='old',action='read')
4372       call getenv_loc('IONPAR',ionname)
4373       open (iion,file=ionname,status='old',action='read')
4374
4375 !      print *,"Processor",myrank," opened file ISIDEP"
4376 !      print *,"Processor",myrank," opened parameter files"
4377 #elif (defined G77)
4378       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4379       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4380 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4381 ! Get parameter filenames and open the parameter files.
4382       call getenv_loc('BONDPAR',bondname)
4383       open (ibond,file=bondname,status='old')
4384       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4385       open (ibond_nucl,file=bondname_nucl,status='old')
4386
4387       call getenv_loc('THETPAR',thetname)
4388       open (ithep,file=thetname,status='old')
4389       call getenv_loc('ROTPAR',rotname)
4390       open (irotam,file=rotname,status='old')
4391       call getenv_loc('TORPAR',torname)
4392       open (itorp,file=torname,status='old')
4393       call getenv_loc('TORDPAR',tordname)
4394       open (itordp,file=tordname,status='old')
4395       call getenv_loc('SCCORPAR',sccorname)
4396       open (isccor,file=sccorname,status='old')
4397       call getenv_loc('FOURIER',fouriername)
4398       open (ifourier,file=fouriername,status='old')
4399       call getenv_loc('ELEPAR',elename)
4400       open (ielep,file=elename,status='old')
4401       call getenv_loc('SIDEPAR',sidename)
4402       open (isidep,file=sidename,status='old')
4403
4404       open (ithep_nucl,file=thetname_nucl,status='old')
4405       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4406       open (irotam_nucl,file=rotname_nucl,status='old')
4407       call getenv_loc('TORPAR_NUCL',torname_nucl)
4408       open (itorp_nucl,file=torname_nucl,status='old')
4409       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4410       open (itordp_nucl,file=tordname_nucl,status='old')
4411       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4412       open (isidep_nucl,file=sidename_nucl,status='old')
4413
4414       call getenv_loc('LIPTRANPAR',liptranname)
4415       open (iliptranpar,file=liptranname,status='old')
4416       call getenv_loc('TUBEPAR',tubename)
4417       open (itube,file=tubename,status='old')
4418       call getenv_loc('IONPAR',ionname)
4419       open (iion,file=ionname,status='old')
4420 #else
4421       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4422         readonly)
4423        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4424 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4425 ! Get parameter filenames and open the parameter files.
4426       call getenv_loc('BONDPAR',bondname)
4427       open (ibond,file=bondname,status='old',action='read')
4428       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4429       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4430       call getenv_loc('THETPAR',thetname)
4431       open (ithep,file=thetname,status='old',action='read')
4432       call getenv_loc('ROTPAR',rotname)
4433       open (irotam,file=rotname,status='old',action='read')
4434       call getenv_loc('TORPAR',torname)
4435       open (itorp,file=torname,status='old',action='read')
4436       call getenv_loc('TORDPAR',tordname)
4437       open (itordp,file=tordname,status='old',action='read')
4438       call getenv_loc('SCCORPAR',sccorname)
4439       open (isccor,file=sccorname,status='old',action='read')
4440 #ifndef CRYST_THETA
4441       call getenv_loc('THETPARPDB',thetname_pdb)
4442       print *,"thetname_pdb ",thetname_pdb
4443       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4444       print *,ithep_pdb," opened"
4445 #endif
4446       call getenv_loc('FOURIER',fouriername)
4447       open (ifourier,file=fouriername,status='old',readonly)
4448       call getenv_loc('ELEPAR',elename)
4449       open (ielep,file=elename,status='old',readonly)
4450       call getenv_loc('SIDEPAR',sidename)
4451       open (isidep,file=sidename,status='old',readonly)
4452
4453       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4454       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4455       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4456       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4457       call getenv_loc('TORPAR_NUCL',torname_nucl)
4458       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4459       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4460       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4461       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4462       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4463       call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4464       open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4465
4466       call getenv_loc('LIPTRANPAR',liptranname)
4467       open (iliptranpar,file=liptranname,status='old',action='read')
4468       call getenv_loc('TUBEPAR',tubename)
4469       open (itube,file=tubename,status='old',action='read')
4470       call getenv_loc('IONPAR',ionname)
4471       open (iion,file=ionname,status='old',action='read')
4472
4473 #ifndef CRYST_SC
4474       call getenv_loc('ROTPARPDB',rotname_pdb)
4475       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4476 #endif
4477 #endif
4478       call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4479 #if defined(WINIFL) || defined(WINPGI)
4480       open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4481 #elif (defined CRAY)  || (defined AIX)
4482       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4483 #elif (defined G77)
4484       open (iscpp_nucl,file=scpname_nucl,status='old')
4485 #else
4486       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4487 #endif
4488
4489 #ifndef OLDSCP
4490 !
4491 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4492 ! Use -DOLDSCP to use hard-coded constants instead.
4493 !
4494       call getenv_loc('SCPPAR',scpname)
4495 #if defined(WINIFL) || defined(WINPGI)
4496       open (iscpp,file=scpname,status='old',readonly,shared)
4497 #elif (defined CRAY)  || (defined AIX)
4498       open (iscpp,file=scpname,status='old',action='read')
4499 #elif (defined G77)
4500       open (iscpp,file=scpname,status='old')
4501 #else
4502       open (iscpp,file=scpname,status='old',action='read')
4503 #endif
4504 #endif
4505       call getenv_loc('PATTERN',patname)
4506 #if defined(WINIFL) || defined(WINPGI)
4507       open (icbase,file=patname,status='old',readonly,shared)
4508 #elif (defined CRAY)  || (defined AIX)
4509       open (icbase,file=patname,status='old',action='read')
4510 #elif (defined G77)
4511       open (icbase,file=patname,status='old')
4512 #else
4513       open (icbase,file=patname,status='old',action='read')
4514 #endif
4515 #ifdef MPI
4516 ! Open output file only for CG processes
4517 !      print *,"Processor",myrank," fg_rank",fg_rank
4518       if (fg_rank.eq.0) then
4519
4520       if (nodes.eq.1) then
4521         npos=3
4522       else
4523         npos = dlog10(dfloat(nodes-1))+1
4524       endif
4525       if (npos.lt.3) npos=3
4526       write (liczba,'(i1)') npos
4527       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4528         //')'
4529       write (liczba,form) me
4530       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4531         liczba(:ilen(liczba))
4532       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4533         //'.int'
4534       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4535         //'.pdb'
4536       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4537         liczba(:ilen(liczba))//'.mol2'
4538       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4539         liczba(:ilen(liczba))//'.stat'
4540       if (lentmp.gt.0) &
4541         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4542             //liczba(:ilen(liczba))//'.stat')
4543       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4544         //'.rst'
4545       if(usampl) then
4546           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4547        liczba(:ilen(liczba))//'.const'
4548       endif
4549
4550       endif
4551 #else
4552       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4553       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4554       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4555       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4556       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4557       if (lentmp.gt.0) &
4558         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4559           '.stat')
4560       rest2name=prefix(:ilen(prefix))//'.rst'
4561       if(usampl) then
4562          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4563       endif
4564 #endif
4565 #if defined(AIX) || defined(PGI)
4566       if (me.eq.king .or. .not. out1file) &
4567          open(iout,file=outname,status='unknown')
4568 #ifdef DEBUG
4569       if (fg_rank.gt.0) then
4570         write (liczba,'(i3.3)') myrank/nfgtasks
4571         write (ll,'(bz,i3.3)') fg_rank
4572         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4573          status='unknown')
4574       endif
4575 #endif
4576       if(me.eq.king) then
4577        open(igeom,file=intname,status='unknown',position='append')
4578        open(ipdb,file=pdbname,status='unknown')
4579        open(imol2,file=mol2name,status='unknown')
4580        open(istat,file=statname,status='unknown',position='append')
4581       else
4582 !1out       open(iout,file=outname,status='unknown')
4583       endif
4584 #else
4585       if (me.eq.king .or. .not.out1file) &
4586           open(iout,file=outname,status='unknown')
4587 #ifdef DEBUG
4588       if (fg_rank.gt.0) then
4589         write (liczba,'(i3.3)') myrank/nfgtasks
4590         write (ll,'(bz,i3.3)') fg_rank
4591         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4592          status='unknown')
4593       endif
4594 #endif
4595       if(me.eq.king) then
4596        open(igeom,file=intname,status='unknown',access='append')
4597        open(ipdb,file=pdbname,status='unknown')
4598        open(imol2,file=mol2name,status='unknown')
4599        open(istat,file=statname,status='unknown',access='append')
4600       else
4601 !1out       open(iout,file=outname,status='unknown')
4602       endif
4603 #endif
4604       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4605       csa_seed=prefix(:lenpre)//'.CSA.seed'
4606       csa_history=prefix(:lenpre)//'.CSA.history'
4607       csa_bank=prefix(:lenpre)//'.CSA.bank'
4608       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4609       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4610       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4611 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4612       csa_int=prefix(:lenpre)//'.int'
4613       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4614       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4615       csa_in=prefix(:lenpre)//'.CSA.in'
4616 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4617 ! Write file names
4618       if (me.eq.king)then
4619       write (iout,'(80(1h-))')
4620       write (iout,'(30x,a)') "FILE ASSIGNMENT"
4621       write (iout,'(80(1h-))')
4622       write (iout,*) "Input file                      : ",&
4623         pref_orig(:ilen(pref_orig))//'.inp'
4624       write (iout,*) "Output file                     : ",&
4625         outname(:ilen(outname))
4626       write (iout,*)
4627       write (iout,*) "Sidechain potential file        : ",&
4628         sidename(:ilen(sidename))
4629 #ifndef OLDSCP
4630       write (iout,*) "SCp potential file              : ",&
4631         scpname(:ilen(scpname))
4632 #endif
4633       write (iout,*) "Electrostatic potential file    : ",&
4634         elename(:ilen(elename))
4635       write (iout,*) "Cumulant coefficient file       : ",&
4636         fouriername(:ilen(fouriername))
4637       write (iout,*) "Torsional parameter file        : ",&
4638         torname(:ilen(torname))
4639       write (iout,*) "Double torsional parameter file : ",&
4640         tordname(:ilen(tordname))
4641       write (iout,*) "SCCOR parameter file : ",&
4642         sccorname(:ilen(sccorname))
4643       write (iout,*) "Bond & inertia constant file    : ",&
4644         bondname(:ilen(bondname))
4645       write (iout,*) "Bending parameter file          : ",&
4646         thetname(:ilen(thetname))
4647       write (iout,*) "Rotamer parameter file          : ",&
4648         rotname(:ilen(rotname))
4649 !el----
4650 #ifndef CRYST_THETA
4651       write (iout,*) "Thetpdb parameter file          : ",&
4652         thetname_pdb(:ilen(thetname_pdb))
4653 #endif
4654 !el
4655       write (iout,*) "Threading database              : ",&
4656         patname(:ilen(patname))
4657       if (lentmp.ne.0) &
4658       write (iout,*)" DIRTMP                          : ",&
4659         tmpdir(:lentmp)
4660       write (iout,'(80(1h-))')
4661       endif
4662       return
4663       end subroutine openunits
4664 !-----------------------------------------------------------------------------
4665       subroutine readrst
4666
4667       use geometry_data, only: nres,dc
4668       use MD_data
4669 !      implicit real*8 (a-h,o-z)
4670 !      include 'DIMENSIONS'
4671 !      include 'COMMON.CHAIN'
4672 !      include 'COMMON.IOUNITS'
4673 !      include 'COMMON.MD'
4674 !el local variables
4675       integer ::i,j
4676 !     real(kind=8) :: var,ene
4677
4678       open(irest2,file=rest2name,status='unknown')
4679       read(irest2,*) totT,EK,potE,totE,t_bath
4680       totTafm=totT
4681 !      do i=1,2*nres
4682 ! AL 4/17/17: Now reading d_t(0,:) too
4683       do i=0,2*nres
4684          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4685       enddo
4686 !      do i=1,2*nres
4687 ! AL 4/17/17: Now reading d_c(0,:) too
4688       do i=0,2*nres
4689          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4690       enddo
4691       if(usampl) then
4692              read (irest2,*) iset
4693       endif
4694       close(irest2)
4695       return
4696       end subroutine readrst
4697 !-----------------------------------------------------------------------------
4698       subroutine read_fragments
4699
4700       use energy_data
4701 !      use geometry
4702       use control_data, only:out1file
4703       use MD_data
4704       use MPI_data
4705 !      implicit real*8 (a-h,o-z)
4706 !      include 'DIMENSIONS'
4707 #ifdef MPI
4708       include 'mpif.h'
4709 #endif
4710 !      include 'COMMON.SETUP'
4711 !      include 'COMMON.CHAIN'
4712 !      include 'COMMON.IOUNITS'
4713 !      include 'COMMON.MD'
4714 !      include 'COMMON.CONTROL'
4715 !el local variables
4716       integer :: i
4717 !     real(kind=8) :: var,ene
4718
4719       read(inp,*) nset,nfrag,npair,nfrag_back
4720
4721 !el from module energy
4722 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
4723       if(.not.allocated(wfrag_back)) then
4724         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4725         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4726
4727         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4728         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
4729
4730         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4731         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
4732       endif
4733
4734       if(me.eq.king.or..not.out1file) &
4735        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4736         " nfrag_back",nfrag_back
4737       do iset=1,nset
4738          read(inp,*) mset(iset)
4739        do i=1,nfrag
4740          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4741            qinfrag(i,iset)
4742          if(me.eq.king.or..not.out1file) &
4743           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4744            ifrag(2,i,iset), qinfrag(i,iset)
4745        enddo
4746        do i=1,npair
4747         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4748           qinpair(i,iset)
4749         if(me.eq.king.or..not.out1file) &
4750          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4751           ipair(2,i,iset), qinpair(i,iset)
4752        enddo
4753        do i=1,nfrag_back
4754         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4755            wfrag_back(3,i,iset),&
4756            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4757         if(me.eq.king.or..not.out1file) &
4758          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4759          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4760        enddo
4761       enddo
4762       return
4763       end subroutine read_fragments
4764 !-----------------------------------------------------------------------------
4765 ! shift.F       io_csa
4766 !-----------------------------------------------------------------------------
4767       subroutine csa_read
4768
4769       use csa_data
4770 !      implicit real*8 (a-h,o-z)
4771 !      include 'DIMENSIONS'
4772 !      include 'COMMON.CSA'
4773 !      include 'COMMON.BANK'
4774 !      include 'COMMON.IOUNITS'
4775 !el local variables
4776 !     integer :: ntf,ik,iw_pdb
4777 !     real(kind=8) :: var,ene
4778
4779       open(icsa_in,file=csa_in,status="old",err=100)
4780        read(icsa_in,*) nconf
4781        read(icsa_in,*) jstart,jend
4782        read(icsa_in,*) nstmax
4783        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4784        read(icsa_in,*) nran0,nran1,irr
4785        read(icsa_in,*) nseed
4786        read(icsa_in,*) ntotal,cut1,cut2
4787        read(icsa_in,*) estop
4788        read(icsa_in,*) icmax,irestart
4789        read(icsa_in,*) ntbankm,dele,difcut
4790        read(icsa_in,*) iref,rmscut,pnccut
4791        read(icsa_in,*) ndiff
4792       close(icsa_in)
4793
4794       return
4795
4796  100  continue
4797       return
4798       end subroutine csa_read
4799 !-----------------------------------------------------------------------------
4800       subroutine initial_write
4801
4802       use csa_data
4803 !      implicit real*8 (a-h,o-z)
4804 !      include 'DIMENSIONS'
4805 !      include 'COMMON.CSA'
4806 !      include 'COMMON.BANK'
4807 !      include 'COMMON.IOUNITS'
4808 !el local variables
4809 !     integer :: ntf,ik,iw_pdb
4810 !     real(kind=8) :: var,ene
4811
4812       open(icsa_seed,file=csa_seed,status="unknown")
4813        write(icsa_seed,*) "seed"
4814       close(31)
4815 #if defined(AIX) || defined(PGI)
4816        open(icsa_history,file=csa_history,status="unknown",&
4817         position="append")
4818 #else
4819        open(icsa_history,file=csa_history,status="unknown",&
4820         access="append")
4821 #endif
4822        write(icsa_history,*) nconf
4823        write(icsa_history,*) jstart,jend
4824        write(icsa_history,*) nstmax
4825        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4826        write(icsa_history,*) nran0,nran1,irr
4827        write(icsa_history,*) nseed
4828        write(icsa_history,*) ntotal,cut1,cut2
4829        write(icsa_history,*) estop
4830        write(icsa_history,*) icmax,irestart
4831        write(icsa_history,*) ntbankm,dele,difcut
4832        write(icsa_history,*) iref,rmscut,pnccut
4833        write(icsa_history,*) ndiff
4834
4835        write(icsa_history,*)
4836       close(icsa_history)
4837
4838       open(icsa_bank1,file=csa_bank1,status="unknown")
4839        write(icsa_bank1,*) 0
4840       close(icsa_bank1)
4841
4842       return
4843       end subroutine initial_write
4844 !-----------------------------------------------------------------------------
4845       subroutine restart_write
4846
4847       use csa_data
4848 !      implicit real*8 (a-h,o-z)
4849 !      include 'DIMENSIONS'
4850 !      include 'COMMON.IOUNITS'
4851 !      include 'COMMON.CSA'
4852 !      include 'COMMON.BANK'
4853 !el local variables
4854 !     integer :: ntf,ik,iw_pdb
4855 !     real(kind=8) :: var,ene
4856
4857 #if defined(AIX) || defined(PGI)
4858        open(icsa_history,file=csa_history,position="append")
4859 #else
4860        open(icsa_history,file=csa_history,access="append")
4861 #endif
4862        write(icsa_history,*)
4863        write(icsa_history,*) "This is restart"
4864        write(icsa_history,*)
4865        write(icsa_history,*) nconf
4866        write(icsa_history,*) jstart,jend
4867        write(icsa_history,*) nstmax
4868        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4869        write(icsa_history,*) nran0,nran1,irr
4870        write(icsa_history,*) nseed
4871        write(icsa_history,*) ntotal,cut1,cut2
4872        write(icsa_history,*) estop
4873        write(icsa_history,*) icmax,irestart
4874        write(icsa_history,*) ntbankm,dele,difcut
4875        write(icsa_history,*) iref,rmscut,pnccut
4876        write(icsa_history,*) ndiff
4877        write(icsa_history,*)
4878        write(icsa_history,*) "irestart is: ", irestart
4879
4880        write(icsa_history,*)
4881       close(icsa_history)
4882
4883       return
4884       end subroutine restart_write
4885 !-----------------------------------------------------------------------------
4886 ! test.F
4887 !-----------------------------------------------------------------------------
4888       subroutine write_pdb(npdb,titelloc,ee)
4889
4890 !      implicit real*8 (a-h,o-z)
4891 !      include 'DIMENSIONS'
4892 !      include 'COMMON.IOUNITS'
4893       character(len=50) :: titelloc1
4894       character*(*) :: titelloc
4895       character(len=3) :: zahl
4896       character(len=5) :: liczba5
4897       real(kind=8) :: ee
4898       integer :: npdb   !,ilen
4899 !el      external ilen
4900 !el local variables
4901       integer :: lenpre
4902 !     real(kind=8) :: var,ene
4903
4904       titelloc1=titelloc
4905       lenpre=ilen(prefix)
4906       if (npdb.lt.1000) then
4907        call numstr(npdb,zahl)
4908        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4909       else
4910         if (npdb.lt.10000) then
4911          write(liczba5,'(i1,i4)') 0,npdb
4912         else
4913          write(liczba5,'(i5)') npdb
4914         endif
4915         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4916       endif
4917       call pdbout(ee,titelloc1,ipdb)
4918       close(ipdb)
4919       return
4920       end subroutine write_pdb
4921 !-----------------------------------------------------------------------------
4922 ! thread.F
4923 !-----------------------------------------------------------------------------
4924       subroutine write_thread_summary
4925 ! Thread the sequence through a database of known structures
4926       use control_data, only: refstr
4927 !      use geometry
4928       use energy_data, only: n_ene_comp
4929       use compare_data
4930 !      implicit real*8 (a-h,o-z)
4931 !      include 'DIMENSIONS'
4932 #ifdef MPI
4933       use MPI_data      !include 'COMMON.INFO'
4934       include 'mpif.h'
4935 #endif
4936 !      include 'COMMON.CONTROL'
4937 !      include 'COMMON.CHAIN'
4938 !      include 'COMMON.DBASE'
4939 !      include 'COMMON.INTERACT'
4940 !      include 'COMMON.VAR'
4941 !      include 'COMMON.THREAD'
4942 !      include 'COMMON.FFIELD'
4943 !      include 'COMMON.SBRIDGE'
4944 !      include 'COMMON.HEADER'
4945 !      include 'COMMON.NAMES'
4946 !      include 'COMMON.IOUNITS'
4947 !      include 'COMMON.TIME1'
4948
4949       integer,dimension(maxthread) :: ip
4950       real(kind=8),dimension(0:n_ene) :: energia
4951 !el local variables
4952       integer :: i,j,ii,jj,ipj,ik,kk,ist
4953       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4954
4955       write (iout,'(30x,a/)') &
4956        '  *********** Summary threading statistics ************'
4957       write (iout,'(a)') 'Initial energies:'
4958       write (iout,'(a4,2x,a12,14a14,3a8)') &
4959        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4960        'RMSnat','NatCONT','NNCONT','RMS'
4961 ! Energy sort patterns
4962       do i=1,nthread
4963         ip(i)=i
4964       enddo
4965       do i=1,nthread-1
4966         enet=ener(n_ene-1,ip(i))
4967         jj=i
4968         do j=i+1,nthread
4969           if (ener(n_ene-1,ip(j)).lt.enet) then
4970             jj=j
4971             enet=ener(n_ene-1,ip(j))
4972           endif
4973         enddo
4974         if (jj.ne.i) then
4975           ipj=ip(jj)
4976           ip(jj)=ip(i)
4977           ip(i)=ipj
4978         endif
4979       enddo
4980       do ik=1,nthread
4981         i=ip(ik)
4982         ii=ipatt(1,i)
4983         ist=nres_base(2,ii)+ipatt(2,i)
4984         do kk=1,n_ene_comp
4985           energia(i)=ener0(kk,i)
4986         enddo
4987         etot=ener0(n_ene_comp+1,i)
4988         rmsnat=ener0(n_ene_comp+2,i)
4989         rms=ener0(n_ene_comp+3,i)
4990         frac=ener0(n_ene_comp+4,i)
4991         frac_nn=ener0(n_ene_comp+5,i)
4992
4993         if (refstr) then
4994         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4995         i,str_nam(ii),ist+1,&
4996         (energia(print_order(kk)),kk=1,nprint_ene),&
4997         etot,rmsnat,frac,frac_nn,rms
4998         else
4999         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5000         i,str_nam(ii),ist+1,&
5001         (energia(print_order(kk)),kk=1,nprint_ene),etot
5002         endif
5003       enddo
5004       write (iout,'(//a)') 'Final energies:'
5005       write (iout,'(a4,2x,a12,17a14,3a8)') &
5006        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5007        'RMSnat','NatCONT','NNCONT','RMS'
5008       do ik=1,nthread
5009         i=ip(ik)
5010         ii=ipatt(1,i)
5011         ist=nres_base(2,ii)+ipatt(2,i)
5012         do kk=1,n_ene_comp
5013           energia(kk)=ener(kk,ik)
5014         enddo
5015         etot=ener(n_ene_comp+1,i)
5016         rmsnat=ener(n_ene_comp+2,i)
5017         rms=ener(n_ene_comp+3,i)
5018         frac=ener(n_ene_comp+4,i)
5019         frac_nn=ener(n_ene_comp+5,i)
5020         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5021         i,str_nam(ii),ist+1,&
5022         (energia(print_order(kk)),kk=1,nprint_ene),&
5023         etot,rmsnat,frac,frac_nn,rms
5024       enddo
5025       write (iout,'(/a/)') 'IEXAM array:'
5026       write (iout,'(i5)') nexcl
5027       do i=1,nexcl
5028         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5029       enddo
5030       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5031        'Max. time for threading step ',max_time_for_thread,&
5032        'Average time for threading step: ',ave_time_for_thread
5033       return
5034       end subroutine write_thread_summary
5035 !-----------------------------------------------------------------------------
5036       subroutine write_stat_thread(ithread,ipattern,ist)
5037
5038       use energy_data, only: n_ene_comp
5039       use compare_data
5040 !      implicit real*8 (a-h,o-z)
5041 !      include "DIMENSIONS"
5042 !      include "COMMON.CONTROL"
5043 !      include "COMMON.IOUNITS"
5044 !      include "COMMON.THREAD"
5045 !      include "COMMON.FFIELD"
5046 !      include "COMMON.DBASE"
5047 !      include "COMMON.NAMES"
5048       real(kind=8),dimension(0:n_ene) :: energia
5049 !el local variables
5050       integer :: ithread,ipattern,ist,i
5051       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5052
5053 #if defined(AIX) || defined(PGI)
5054       open(istat,file=statname,position='append')
5055 #else
5056       open(istat,file=statname,access='append')
5057 #endif
5058       do i=1,n_ene_comp
5059         energia(i)=ener(i,ithread)
5060       enddo
5061       etot=ener(n_ene_comp+1,ithread)
5062       rmsnat=ener(n_ene_comp+2,ithread)
5063       rms=ener(n_ene_comp+3,ithread)
5064       frac=ener(n_ene_comp+4,ithread)
5065       frac_nn=ener(n_ene_comp+5,ithread)
5066       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5067         ithread,str_nam(ipattern),ist+1,&
5068         (energia(print_order(i)),i=1,nprint_ene),&
5069         etot,rmsnat,frac,frac_nn,rms
5070       close (istat)
5071       return
5072       end subroutine write_stat_thread
5073 !-----------------------------------------------------------------------------
5074 #endif
5075 !-----------------------------------------------------------------------------
5076       end module io_config