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
2407       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2408       bad(:,:)=0.0D0
2409
2410 #ifdef OLDSCP
2411 !
2412 ! Define the SC-p interaction constants (hard-coded; old style)
2413 !
2414       do i=1,ntyp
2415 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2416 ! helix formation)
2417 !       aad(i,1)=0.3D0*4.0D0**12
2418 ! Following line for constants currently implemented
2419 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2420         aad(i,1)=1.5D0*4.0D0**12
2421 !       aad(i,1)=0.17D0*5.6D0**12
2422         aad(i,2)=aad(i,1)
2423 ! "Soft" SC-p repulsion
2424         bad(i,1)=0.0D0
2425 ! Following line for constants currently implemented
2426 !       aad(i,1)=0.3D0*4.0D0**6
2427 ! "Hard" SC-p repulsion
2428         bad(i,1)=3.0D0*4.0D0**6
2429 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2430         bad(i,2)=bad(i,1)
2431 !       aad(i,1)=0.0D0
2432 !       aad(i,2)=0.0D0
2433 !       bad(i,1)=1228.8D0
2434 !       bad(i,2)=1228.8D0
2435       enddo
2436 #else
2437 !
2438 ! 8/9/01 Read the SC-p interaction constants from file
2439 !
2440       do i=1,ntyp
2441         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2442       enddo
2443       do i=1,ntyp
2444         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2445         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2446         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2447         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2448       enddo
2449 !      lprint=.true.
2450       if (lprint) then
2451         write (iout,*) "Parameters of SC-p interactions:"
2452         do i=1,ntyp
2453           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2454            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2455         enddo
2456       endif
2457 !      lprint=.false.
2458 #endif
2459       allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2460
2461       do i=1,ntyp_molec(2)
2462         read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2463       enddo
2464       do i=1,ntyp_molec(2)
2465         aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2466         bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2467       enddo
2468       r0pp=1.12246204830937298142*5.16158
2469       epspp=4.95713/4
2470       AEES=108.661
2471       BEES=0.433246
2472
2473 !
2474 ! Define the constants of the disulfide bridge
2475 !
2476       ebr=-5.50D0
2477 !
2478 ! Old arbitrary potential - commented out.
2479 !
2480 !      dbr= 4.20D0
2481 !      fbr= 3.30D0
2482 !
2483 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2484 ! energy surface of diethyl disulfide.
2485 ! A. Liwo and U. Kozlowska, 11/24/03
2486 !
2487       D0CM = 3.78d0
2488       AKCM = 15.1d0
2489       AKTH = 11.0d0
2490       AKCT = 12.0d0
2491       V1SS =-1.08d0
2492       V2SS = 7.61d0
2493       V3SS = 13.7d0
2494 !      akcm=0.0d0
2495 !      akth=0.0d0
2496 !      akct=0.0d0
2497 !      v1ss=0.0d0
2498 !      v2ss=0.0d0
2499 !      v3ss=0.0d0
2500
2501       if(me.eq.king) then
2502       write (iout,'(/a)') "Disulfide bridge parameters:"
2503       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2504       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2505       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2506       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2507         ' v3ss:',v3ss
2508       endif
2509       return
2510   111 write (iout,*) "Error reading bending energy parameters."
2511       goto 999
2512   112 write (iout,*) "Error reading rotamer energy parameters."
2513       goto 999
2514   113 write (iout,*) "Error reading torsional energy parameters."
2515       goto 999
2516   114 write (iout,*) "Error reading double torsional energy parameters."
2517       goto 999
2518   115 write (iout,*) &
2519         "Error reading cumulant (multibody energy) parameters."
2520       goto 999
2521   116 write (iout,*) "Error reading electrostatic energy parameters."
2522       goto 999
2523   117 write (iout,*) "Error reading side chain interaction parameters."
2524       goto 999
2525   118 write (iout,*) "Error reading SCp interaction parameters."
2526       goto 999
2527   119 write (iout,*) "Error reading SCCOR parameters"
2528   999 continue
2529 #ifdef MPI
2530       call MPI_Finalize(Ierror)
2531 #endif
2532       stop
2533       return
2534       end subroutine parmread
2535 #endif
2536 !-----------------------------------------------------------------------------
2537 ! printmat.f
2538 !-----------------------------------------------------------------------------
2539       subroutine printmat(ldim,m,n,iout,key,a)
2540
2541       integer :: n,ldim
2542       character(len=3),dimension(n) :: key
2543       real(kind=8),dimension(ldim,n) :: a
2544 !el local variables
2545       integer :: i,j,k,m,iout,nlim
2546
2547       do 1 i=1,n,8
2548       nlim=min0(i+7,n)
2549       write (iout,1000) (key(k),k=i,nlim)
2550       write (iout,1020)
2551  1000 format (/5x,8(6x,a3))
2552  1020 format (/80(1h-)/)
2553       do 2 j=1,n
2554       write (iout,1010) key(j),(a(j,k),k=i,nlim)
2555     2 continue
2556     1 continue
2557  1010 format (a3,2x,8(f9.4))
2558       return
2559       end subroutine printmat
2560 !-----------------------------------------------------------------------------
2561 ! readpdb.F
2562 !-----------------------------------------------------------------------------
2563       subroutine readpdb
2564 ! Read the PDB file and convert the peptide geometry into virtual-chain
2565 ! geometry.
2566       use geometry_data
2567       use energy_data, only: itype,istype
2568       use control_data
2569       use compare_data
2570       use MPI_data
2571       use control, only: rescode,sugarcode
2572 !      implicit real*8 (a-h,o-z)
2573 !      include 'DIMENSIONS'
2574 !      include 'COMMON.LOCAL'
2575 !      include 'COMMON.VAR'
2576 !      include 'COMMON.CHAIN'
2577 !      include 'COMMON.INTERACT'
2578 !      include 'COMMON.IOUNITS'
2579 !      include 'COMMON.GEO'
2580 !      include 'COMMON.NAMES'
2581 !      include 'COMMON.CONTROL'
2582 !      include 'COMMON.DISTFIT'
2583 !      include 'COMMON.SETUP'
2584       integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2585 !        ishift_pdb
2586       logical :: lprn=.true.,fail
2587       real(kind=8),dimension(3) :: e1,e2,e3
2588       real(kind=8) :: dcj,efree_temp
2589       character(len=3) :: seq,res
2590       character(len=5) :: atom
2591       character(len=80) :: card
2592       real(kind=8),dimension(3,20) :: sccor
2593       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin       !rescode,
2594       integer :: isugar,molecprev,firstion
2595       character*1 :: sugar
2596       real(kind=8) :: cou
2597       real(kind=8),dimension(3) :: ccc
2598 !el local varables
2599       integer,dimension(2,maxres/3) :: hfrag_alloc
2600       integer,dimension(4,maxres/3) :: bfrag_alloc
2601       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2602       real(kind=8),dimension(:,:), allocatable  :: c_temporary
2603       integer,dimension(:,:) , allocatable  :: itype_temporary
2604       integer,dimension(:),allocatable :: istype_temp
2605       efree_temp=0.0d0
2606       ibeg=1
2607       ishift1=0
2608       ishift=0
2609       molecule=1
2610       counter=0
2611 !      write (2,*) "UNRES_PDB",unres_pdb
2612       ires=0
2613       ires_old=0
2614       nres=0
2615       iii=0
2616       lsecondary=.false.
2617       nhfrag=0
2618       nbfrag=0
2619 !-----------------------------
2620       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2621       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2622
2623       do i=1,100000
2624         read (ipdbin,'(a80)',end=10) card
2625 !       write (iout,'(a)') card
2626         if (card(:5).eq.'HELIX') then
2627           nhfrag=nhfrag+1
2628           lsecondary=.true.
2629           read(card(22:25),*) hfrag(1,nhfrag)
2630           read(card(34:37),*) hfrag(2,nhfrag)
2631         endif
2632         if (card(:5).eq.'SHEET') then
2633           nbfrag=nbfrag+1
2634           lsecondary=.true.
2635           read(card(24:26),*) bfrag(1,nbfrag)
2636           read(card(35:37),*) bfrag(2,nbfrag)
2637 !rc----------------------------------------
2638 !rc  to be corrected !!!
2639           bfrag(3,nbfrag)=bfrag(1,nbfrag)
2640           bfrag(4,nbfrag)=bfrag(2,nbfrag)
2641 !rc----------------------------------------
2642         endif
2643         if (card(:3).eq.'END') then
2644           goto 10
2645         else if (card(:3).eq.'TER') then
2646 ! End current chain
2647           ires_old=ires+2
2648           ishift=ishift+1
2649           ishift1=ishift1+1
2650           itype(ires_old,molecule)=ntyp1_molec(molecule)
2651           itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2652           nres_molec(molecule)=nres_molec(molecule)+2
2653           ibeg=2
2654 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2655           if (unres_pdb) then
2656             do j=1,3
2657               dc(j,ires)=sccor(j,iii)
2658             enddo
2659           else
2660             call sccenter(ires,iii,sccor)
2661 !          iii=0
2662           endif
2663           iii=0
2664         endif
2665 ! Read free energy
2666         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2667 ! Fish out the ATOM cards.
2668 !        write(iout,*) 'card',card(1:20)
2669         if (index(card(1:4),'ATOM').gt.0) then
2670           read (card(12:16),*) atom
2671 !          write (iout,*) "! ",atom," !",ires
2672 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2673           read (card(23:26),*) ires
2674           read (card(18:20),'(a3)') res
2675 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2676 !     &      " ires_old",ires_old
2677 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2678 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2679           if (ires-ishift+ishift1.ne.ires_old) then
2680 ! Calculate the CM of the preceding residue.
2681 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2682             if (ibeg.eq.0) then
2683 !              write (iout,*) "Calculating sidechain center iii",iii
2684               if (unres_pdb) then
2685                 do j=1,3
2686                   dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2687                 enddo
2688               else
2689                 call sccenter(ires_old,iii,sccor)
2690               endif
2691               iii=0
2692             endif
2693 ! Start new residue.
2694             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2695               ires=ires_old
2696               cycle
2697             else if (ibeg.eq.1) then
2698               write (iout,*) "BEG ires",ires
2699               ishift=ires-1
2700               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2701                 ishift=ishift-1
2702                 itype(1,1)=ntyp1
2703                 nres_molec(molecule)=nres_molec(molecule)+1
2704               endif
2705               ires=ires-ishift+ishift1
2706               ires_old=ires
2707 !              write (iout,*) "ishift",ishift," ires",ires,&
2708 !               " ires_old",ires_old
2709               ibeg=0
2710             else if (ibeg.eq.2) then
2711 ! Start a new chain
2712               ishift=-ires_old+ires-1 !!!!!
2713               ishift1=ishift1-1    !!!!!
2714 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2715               ires=ires-ishift+ishift1
2716 !              print *,ires,ishift,ishift1
2717               ires_old=ires
2718               ibeg=0
2719             else
2720               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2721               ires=ires-ishift+ishift1
2722               ires_old=ires
2723             endif
2724 !            print *,'atom',ires,atom
2725             if (res.eq.'ACE' .or. res.eq.'NHE') then
2726               itype(ires,1)=10
2727             else
2728              if (atom.eq.'CA  '.or.atom.eq.'N   ') then
2729              molecule=1
2730               itype(ires,molecule)=rescode(ires,res,0,molecule)
2731               firstion=0
2732 !              nres_molec(molecule)=nres_molec(molecule)+1
2733             else
2734              molecule=2
2735               itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2736 !              nres_molec(molecule)=nres_molec(molecule)+1
2737              read (card(19:19),'(a1)') sugar
2738              isugar=sugarcode(sugar,ires)
2739 !            if (ibeg.eq.1) then
2740 !              istype(1)=isugar
2741 !            else
2742               istype(ires)=isugar
2743 !              print *,"ires=",ires,istype(ires)
2744 !            endif
2745
2746             endif
2747             endif
2748           else
2749             ires=ires-ishift+ishift1
2750           endif
2751 !          write (iout,*) "ires_old",ires_old," ires",ires
2752           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2753 !            ishift1=ishift1+1
2754           endif
2755 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2756           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2757              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2758             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2759 !              print *,ires,ishift,ishift1
2760 !            write (iout,*) "backbone ",atom
2761 #ifdef DEBUG
2762             write (iout,'(2i3,2x,a,3f8.3)') &
2763             ires,itype(ires,1),res,(c(j,ires),j=1,3)
2764 #endif
2765             iii=iii+1
2766               nres_molec(molecule)=nres_molec(molecule)+1
2767             do j=1,3
2768               sccor(j,iii)=c(j,ires)
2769             enddo
2770           else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2771                atom.eq."C2'" .or. atom.eq."C3'" &
2772                .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2773             read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2774 !c            write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2775 !              print *,ires,ishift,ishift1
2776             counter=counter+1
2777 !            iii=iii+1
2778 !            do j=1,3
2779 !              sccor(j,iii)=c(j,ires)
2780 !            enddo
2781             do j=1,3
2782               c(j,ires)=c(j,ires)+ccc(j)/5.0
2783             enddo
2784              print *,counter,molecule
2785              if (counter.eq.5) then
2786 !            iii=iii+1
2787               nres_molec(molecule)=nres_molec(molecule)+1
2788               firstion=0
2789 !            do j=1,3
2790 !              sccor(j,iii)=c(j,ires)
2791 !            enddo
2792              counter=0
2793            endif
2794 !            print *, "ATOM",atom(1:3)
2795           else if (atom.eq."C5'") then
2796              read (card(19:19),'(a1)') sugar
2797              isugar=sugarcode(sugar,ires)
2798             if (ibeg.eq.1) then
2799               istype(1)=isugar
2800             else
2801               istype(ires)=isugar
2802 !              print *,ires,istype(ires)
2803             endif
2804             if (unres_pdb) then
2805               read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2806             else
2807               iii=iii+1
2808               read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2809             endif
2810 !            write (*,*) card(23:27),ires,itype(ires,1)
2811           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2812                    atom.ne.'N' .and. atom.ne.'C' .and. &
2813                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2814                    atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2815                    .and. atom.ne.'P  '.and. &
2816                   atom(1:1).ne.'H' .and. &
2817                   atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2818                   ) then
2819 !            write (iout,*) "sidechain ",atom
2820 !            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2821                  if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2822 !                        write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2823             iii=iii+1
2824             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2825               endif
2826           endif
2827         else if ((ions).and.(card(1:6).eq.'HETATM')) then
2828        if (firstion.eq.0) then
2829        firstion=1
2830        if (unres_pdb) then
2831          do j=1,3
2832            dc(j,ires)=sccor(j,iii)
2833          enddo
2834        else
2835           call sccenter(ires,iii,sccor)
2836        endif
2837        endif
2838           read (card(12:16),*) atom
2839           print *,"HETATOM", atom
2840           read (card(18:20),'(a3)') res
2841           if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
2842           (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
2843           .or.(atom(1:2).eq.'K ')) &
2844           then
2845            ires=ires+1
2846            if (molecule.ne.5) molecprev=molecule
2847            molecule=5
2848            nres_molec(molecule)=nres_molec(molecule)+1
2849            print *,"HERE",nres_molec(molecule)
2850            itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2851            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2852           endif
2853         endif !atom
2854       enddo
2855    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2856       if (ires.eq.0) return
2857 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2858 ! system
2859       nres=ires
2860       if ((ires_old.ne.ires).and.(molecule.ne.5)) &
2861          nres_molec(molecule)=nres_molec(molecule)-2
2862 !      print *,'I have', nres_molec(:)
2863
2864       do k=1,4 ! ions are without dummy
2865        if (nres_molec(k).eq.0) cycle
2866       do i=2,nres-1
2867 !        write (iout,*) i,itype(i,1)
2868 !        if (itype(i,1).eq.ntyp1) then
2869 !          write (iout,*) "dummy",i,itype(i,1)
2870 !          do j=1,3
2871 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2872 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
2873 !            dc(j,i)=c(j,i)
2874 !          enddo
2875 !        endif
2876         if (itype(i,k).eq.ntyp1_molec(k)) then
2877          if (itype(i+1,k).eq.ntyp1_molec(k)) then
2878           if (itype(i+2,k).eq.0) then
2879 !           print *,"masz sieczke"
2880            do j=1,5
2881             if (itype(i+2,j).ne.0) then
2882             itype(i+1,k)=0
2883             itype(i+1,j)=ntyp1_molec(j)
2884             nres_molec(k)=nres_molec(k)-1
2885             nres_molec(j)=nres_molec(j)+1
2886             go to 3331
2887             endif
2888            enddo
2889  3331      continue
2890           endif
2891 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2892 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
2893 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
2894            if (unres_pdb) then
2895 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2896 !            print *,i,'tu dochodze'
2897             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2898             if (fail) then
2899               e2(1)=0.0d0
2900               e2(2)=1.0d0
2901               e2(3)=0.0d0
2902             endif !fail
2903 !            print *,i,'a tu?'
2904             do j=1,3
2905              c(j,i)=c(j,i-1)-1.9d0*e2(j)
2906             enddo
2907            else   !unres_pdb
2908            do j=1,3
2909              dcj=(c(j,i-2)-c(j,i-3))/2.0
2910             if (dcj.eq.0) dcj=1.23591524223
2911              c(j,i)=c(j,i-1)+dcj
2912              c(j,nres+i)=c(j,i)
2913            enddo
2914           endif   !unres_pdb
2915          else     !itype(i+1,1).eq.ntyp1
2916           if (unres_pdb) then
2917 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2918             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2919             if (fail) then
2920               e2(1)=0.0d0
2921               e2(2)=1.0d0
2922               e2(3)=0.0d0
2923             endif
2924             do j=1,3
2925               c(j,i)=c(j,i+1)-1.9d0*e2(j)
2926             enddo
2927           else !unres_pdb
2928            do j=1,3
2929             dcj=(c(j,i+3)-c(j,i+2))/2.0
2930             if (dcj.eq.0) dcj=1.23591524223
2931             c(j,i)=c(j,i+1)-dcj
2932             c(j,nres+i)=c(j,i)
2933            enddo
2934           endif !unres_pdb
2935          endif !itype(i+1,1).eq.ntyp1
2936         endif  !itype.eq.ntyp1
2937
2938       enddo
2939       enddo
2940 ! Calculate the CM of the last side chain.
2941       if (iii.gt.0)  then
2942       if (unres_pdb) then
2943         do j=1,3
2944           dc(j,ires)=sccor(j,iii)
2945         enddo
2946       else
2947         call sccenter(ires,iii,sccor)
2948       endif
2949       endif
2950 !      nres=ires
2951       nsup=nres
2952       nstart_sup=1
2953 !      print *,"molecule",molecule
2954       if ((itype(nres,1).ne.10)) then
2955         nres=nres+1
2956           if (molecule.eq.5) molecule=molecprev
2957         itype(nres,molecule)=ntyp1_molec(molecule)
2958         nres_molec(molecule)=nres_molec(molecule)+1
2959         if (unres_pdb) then
2960 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2961           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2962           if (fail) then
2963             e2(1)=0.0d0
2964             e2(2)=1.0d0
2965             e2(3)=0.0d0
2966           endif
2967           do j=1,3
2968             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2969           enddo
2970         else
2971         do j=1,3
2972           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2973           c(j,nres)=c(j,nres-1)+dcj
2974           c(j,2*nres)=c(j,nres)
2975         enddo
2976         endif
2977       endif
2978 !     print *,'I have',nres, nres_molec(:)
2979
2980 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2981 #ifdef WHAM_RUN
2982       if (nres.ne.nres0) then
2983         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2984                        " NRES0=",nres0
2985         stop "Error nres value in WHAM input"
2986       endif
2987 #endif
2988 !---------------------------------
2989 !el reallocate tables
2990 !      do i=1,maxres/3
2991 !       do j=1,2
2992 !         hfrag_alloc(j,i)=hfrag(j,i)
2993 !        enddo
2994 !       do j=1,4
2995 !         bfrag_alloc(j,i)=bfrag(j,i)
2996 !        enddo
2997 !      enddo
2998
2999 !      deallocate(hfrag)
3000 !      deallocate(bfrag)
3001 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
3002 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3003 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3004 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
3005
3006 !      do i=1,nhfrag
3007 !       do j=1,2
3008 !         hfrag(j,i)=hfrag_alloc(j,i)
3009 !        enddo
3010 !      enddo
3011 !      do i=1,nbfrag
3012 !       do j=1,4
3013 !         bfrag(j,i)=bfrag_alloc(j,i)
3014 !        enddo
3015 !      enddo
3016 !el end reallocate tables
3017 !---------------------------------
3018       do i=2,nres-1
3019         do j=1,3
3020           c(j,i+nres)=dc(j,i)
3021         enddo
3022       enddo
3023       do j=1,3
3024         c(j,nres+1)=c(j,1)
3025         c(j,2*nres)=c(j,nres)
3026       enddo
3027
3028       if (itype(1,1).eq.ntyp1) then
3029         nsup=nsup-1
3030         nstart_sup=2
3031         if (unres_pdb) then
3032 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3033           call refsys(2,3,4,e1,e2,e3,fail)
3034           if (fail) then
3035             e2(1)=0.0d0
3036             e2(2)=1.0d0
3037             e2(3)=0.0d0
3038           endif
3039           do j=1,3
3040             c(j,1)=c(j,2)-1.9d0*e2(j)
3041           enddo
3042         else
3043         do j=1,3
3044           dcj=(c(j,4)-c(j,3))/2.0
3045           c(j,1)=c(j,2)-dcj
3046           c(j,nres+1)=c(j,1)
3047         enddo
3048         endif
3049       endif
3050 ! First lets assign correct dummy to correct type of chain
3051 ! 1) First residue
3052       if (itype(1,1).eq.ntyp1) then
3053         if (itype(2,1).eq.0) then
3054          do j=2,5
3055            if (itype(2,j).ne.0) then
3056            itype(1,1)=0
3057            itype(1,j)=ntyp1_molec(j)
3058            nres_molec(1)=nres_molec(1)-1
3059            nres_molec(j)=nres_molec(j)+1
3060            go to 3231
3061            endif
3062          enddo
3063 3231    continue
3064         endif
3065        endif
3066        print *,'I have',nres, nres_molec(:)
3067
3068 ! Copy the coordinates to reference coordinates
3069 !      do i=1,2*nres
3070 !        do j=1,3
3071 !          cref(j,i)=c(j,i)
3072 !        enddo
3073 !      enddo
3074 ! Calculate internal coordinates.
3075       if (lprn) then
3076       write (iout,'(/a)') &
3077         "Cartesian coordinates of the reference structure"
3078       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3079        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3080       do ires=1,nres
3081         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3082           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3083           (c(j,ires+nres),j=1,3)
3084       enddo
3085       endif
3086 ! znamy już nres wiec mozna alokowac tablice
3087 ! Calculate internal coordinates.
3088       if(me.eq.king.or..not.out1file)then
3089        write (iout,'(a)') &
3090          "Backbone and SC coordinates as read from the PDB"
3091        do ires=1,nres
3092         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3093           ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3094           (c(j,nres+ires),j=1,3)
3095        enddo
3096       endif
3097 ! NOW LETS ROCK! SORTING
3098       allocate(c_temporary(3,2*nres))
3099       allocate(itype_temporary(nres,5))
3100       allocate(molnum(nres))
3101       allocate(istype_temp(nres))
3102        itype_temporary(:,:)=0
3103       seqalingbegin=1
3104       do k=1,5
3105         do i=1,nres
3106          if (itype(i,k).ne.0) then
3107           do j=1,3
3108           c_temporary(j,seqalingbegin)=c(j,i)
3109           c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3110
3111           enddo
3112           itype_temporary(seqalingbegin,k)=itype(i,k)
3113           print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3114           istype_temp(seqalingbegin)=istype(i)
3115           molnum(seqalingbegin)=k
3116           seqalingbegin=seqalingbegin+1
3117          endif
3118         enddo
3119        enddo
3120        do i=1,2*nres
3121         do j=1,3
3122         c(j,i)=c_temporary(j,i)
3123         enddo
3124        enddo
3125        do k=1,5
3126         do i=1,nres
3127          itype(i,k)=itype_temporary(i,k)
3128          istype(i)=istype_temp(i)
3129         enddo
3130        enddo
3131       if (itype(1,1).eq.ntyp1) then
3132         nsup=nsup-1
3133         nstart_sup=2
3134         if (unres_pdb) then
3135 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3136           call refsys(2,3,4,e1,e2,e3,fail)
3137           if (fail) then
3138             e2(1)=0.0d0
3139             e2(2)=1.0d0
3140             e2(3)=0.0d0
3141           endif
3142           do j=1,3
3143             c(j,1)=c(j,2)-1.9d0*e2(j)
3144           enddo
3145         else
3146         do j=1,3
3147           dcj=(c(j,4)-c(j,3))/2.0
3148           c(j,1)=c(j,2)-dcj
3149           c(j,nres+1)=c(j,1)
3150         enddo
3151         endif
3152       endif
3153
3154       if (lprn) then
3155       write (iout,'(/a)') &
3156         "Cartesian coordinates of the reference structure after sorting"
3157       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3158        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3159       do ires=1,nres
3160         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3161           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3162           (c(j,ires+nres),j=1,3)
3163       enddo
3164       endif
3165
3166 !       print *,seqalingbegin,nres
3167       if(.not.allocated(vbld)) then
3168        allocate(vbld(2*nres))
3169        do i=1,2*nres
3170          vbld(i)=0.d0
3171        enddo
3172       endif
3173       if(.not.allocated(vbld_inv)) then
3174        allocate(vbld_inv(2*nres))
3175        do i=1,2*nres
3176          vbld_inv(i)=0.d0
3177        enddo
3178       endif
3179 !!!el
3180       if(.not.allocated(theta)) then
3181         allocate(theta(nres+2))
3182         theta(:)=0.0d0
3183       endif
3184
3185       if(.not.allocated(phi)) allocate(phi(nres+2))
3186       if(.not.allocated(alph)) allocate(alph(nres+2))
3187       if(.not.allocated(omeg)) allocate(omeg(nres+2))
3188       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3189       if(.not.allocated(phiref)) allocate(phiref(nres+2))
3190       if(.not.allocated(costtab)) allocate(costtab(nres))
3191       if(.not.allocated(sinttab)) allocate(sinttab(nres))
3192       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3193       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3194       if(.not.allocated(xxref)) allocate(xxref(nres))
3195       if(.not.allocated(yyref)) allocate(yyref(nres))
3196       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3197       if(.not.allocated(dc_norm)) then
3198 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3199         allocate(dc_norm(3,0:2*nres+2))
3200         dc_norm(:,:)=0.d0
3201       endif
3202
3203       call int_from_cart(.true.,.false.)
3204       call sc_loc_geom(.false.)
3205       do i=1,nres
3206         thetaref(i)=theta(i)
3207         phiref(i)=phi(i)
3208       enddo
3209 !      do i=1,2*nres
3210 !        vbld_inv(i)=0.d0
3211 !        vbld(i)=0.d0
3212 !      enddo
3213
3214       do i=1,nres-1
3215         do j=1,3
3216           dc(j,i)=c(j,i+1)-c(j,i)
3217           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3218         enddo
3219       enddo
3220       do i=2,nres-1
3221         do j=1,3
3222           dc(j,i+nres)=c(j,i+nres)-c(j,i)
3223           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3224         enddo
3225 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3226 !        vbld_inv(i+nres)
3227       enddo
3228 !      call chainbuild
3229 ! Copy the coordinates to reference coordinates
3230 ! Splits to single chain if occurs
3231
3232 !      do i=1,2*nres
3233 !        do j=1,3
3234 !          cref(j,i,cou)=c(j,i)
3235 !        enddo
3236 !      enddo
3237 !
3238       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3239       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3240       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3241 !-----------------------------
3242       kkk=1
3243       lll=0
3244       cou=1
3245         write (iout,*) "symetr", symetr
3246       do i=1,nres
3247       lll=lll+1
3248 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3249       if (i.gt.1) then
3250       if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3251       chain_length=lll-1
3252       kkk=kkk+1
3253 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3254       lll=1
3255       endif
3256       endif
3257         do j=1,3
3258           cref(j,i,cou)=c(j,i)
3259           cref(j,i+nres,cou)=c(j,i+nres)
3260           if (i.le.nres) then
3261           chain_rep(j,lll,kkk)=c(j,i)
3262           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3263           endif
3264          enddo
3265       enddo
3266       write (iout,*) chain_length
3267       if (chain_length.eq.0) chain_length=nres
3268       do j=1,3
3269       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3270       chain_rep(j,chain_length+nres,symetr) &
3271       =chain_rep(j,chain_length+nres,1)
3272       enddo
3273 ! diagnostic
3274 !       write (iout,*) "spraw lancuchy",chain_length,symetr
3275 !       do i=1,4
3276 !         do kkk=1,chain_length
3277 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3278 !         enddo
3279 !        enddo
3280 ! enddiagnostic
3281 ! makes copy of chains
3282         write (iout,*) "symetr", symetr
3283       do j=1,3
3284       dc(j,0)=c(j,1)
3285       enddo
3286
3287       if (symetr.gt.1) then
3288        call permut(symetr)
3289        nperm=1
3290        do i=1,symetr
3291        nperm=nperm*i
3292        enddo
3293        do i=1,nperm
3294        write(iout,*) (tabperm(i,kkk),kkk=1,4)
3295        enddo
3296        do i=1,nperm
3297         cou=0
3298         do kkk=1,symetr
3299          icha=tabperm(i,kkk)
3300 !         write (iout,*) i,icha
3301          do lll=1,chain_length
3302           cou=cou+1
3303            if (cou.le.nres) then
3304            do j=1,3
3305             kupa=mod(lll,chain_length)
3306             iprzes=(kkk-1)*chain_length+lll
3307             if (kupa.eq.0) kupa=chain_length
3308 !            write (iout,*) "kupa", kupa
3309             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3310             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3311           enddo
3312           endif
3313          enddo
3314         enddo
3315        enddo
3316        endif
3317 !-koniec robienia kopii
3318 ! diag
3319       do kkk=1,nperm
3320       write (iout,*) "nowa struktura", nperm
3321       do i=1,nres
3322         write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3323       cref(2,i,kkk),&
3324       cref(3,i,kkk),cref(1,nres+i,kkk),&
3325       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3326       enddo
3327   100 format (//'              alpha-carbon coordinates       ',&
3328                 '     centroid coordinates'/ &
3329                 '       ', 6X,'X',11X,'Y',11X,'Z', &
3330                                 10X,'X',11X,'Y',11X,'Z')
3331   110 format (a,'(',i3,')',6f12.5)
3332
3333       enddo
3334 !c enddiag
3335       do j=1,nbfrag
3336         do i=1,4
3337          bfrag(i,j)=bfrag(i,j)-ishift
3338         enddo
3339       enddo
3340
3341       do j=1,nhfrag
3342         do i=1,2
3343          hfrag(i,j)=hfrag(i,j)-ishift
3344         enddo
3345       enddo
3346       ishift_pdb=ishift
3347
3348       return
3349       end subroutine readpdb
3350 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3351 !-----------------------------------------------------------------------------
3352 ! readrtns_CSA.F
3353 !-----------------------------------------------------------------------------
3354       subroutine read_control
3355 !
3356 ! Read contorl data
3357 !
3358 !      use geometry_data
3359       use comm_machsw
3360       use energy_data
3361       use control_data
3362       use compare_data
3363       use MCM_data
3364       use map_data
3365       use csa_data
3366       use MD_data
3367       use MPI_data
3368       use random, only: random_init
3369 !      implicit real*8 (a-h,o-z)
3370 !      include 'DIMENSIONS'
3371 #ifdef MP
3372       use prng, only:prng_restart
3373       include 'mpif.h'
3374       logical :: OKRandom!, prng_restart
3375       real(kind=8) :: r1
3376 #endif
3377 !      include 'COMMON.IOUNITS'
3378 !      include 'COMMON.TIME1'
3379 !      include 'COMMON.THREAD'
3380 !      include 'COMMON.SBRIDGE'
3381 !      include 'COMMON.CONTROL'
3382 !      include 'COMMON.MCM'
3383 !      include 'COMMON.MAP'
3384 !      include 'COMMON.HEADER'
3385 !      include 'COMMON.CSA'
3386 !      include 'COMMON.CHAIN'
3387 !      include 'COMMON.MUCA'
3388 !      include 'COMMON.MD'
3389 !      include 'COMMON.FFIELD'
3390 !      include 'COMMON.INTERACT'
3391 !      include 'COMMON.SETUP'
3392 !el      integer :: KDIAG,ICORFL,IXDR
3393 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3394       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3395         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
3396 !      character(len=80) :: ucase
3397       character(len=640) :: controlcard
3398
3399       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3400       integer i
3401
3402       nglob_csa=0
3403       eglob_csa=1d99
3404       nmin_csa=0
3405       read (INP,'(a)') titel
3406       call card_concat(controlcard,.true.)
3407 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3408 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3409       call reada(controlcard,'SEED',seed,0.0D0)
3410       call random_init(seed)
3411 ! Set up the time limit (caution! The time must be input in minutes!)
3412       read_cart=index(controlcard,'READ_CART').gt.0
3413       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3414       call readi(controlcard,'SYM',symetr,1)
3415       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3416       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3417       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3418       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3419       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3420       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3421       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3422       call reada(controlcard,'DRMS',drms,0.1D0)
3423       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3424        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3425        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3426        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3427        write (iout,'(a,f10.1)')'DRMS    = ',drms
3428        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3429        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3430       endif
3431       call readi(controlcard,'NZ_START',nz_start,0)
3432       call readi(controlcard,'NZ_END',nz_end,0)
3433       call readi(controlcard,'IZ_SC',iz_sc,0)
3434       timlim=60.0D0*timlim
3435       safety = 60.0d0*safety
3436       timem=timlim
3437       modecalc=0
3438       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3439 !C SHIELD keyword sets if the shielding effect of side-chains is used
3440 !C 0 denots no shielding is used all peptide are equally despite the
3441 !C solvent accesible area
3442 !C 1 the newly introduced function
3443 !C 2 reseved for further possible developement
3444       call readi(controlcard,'SHIELD',shield_mode,0)
3445 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3446         write(iout,*) "shield_mode",shield_mode
3447 !C  Varibles set size of box
3448       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3449       protein=index(controlcard,"PROTEIN").gt.0
3450       ions=index(controlcard,"IONS").gt.0
3451       nucleic=index(controlcard,"NUCLEIC").gt.0
3452       write (iout,*) "with_theta_constr ",with_theta_constr
3453       AFMlog=(index(controlcard,'AFM'))
3454       selfguide=(index(controlcard,'SELFGUIDE'))
3455       print *,'AFMlog',AFMlog,selfguide,"KUPA"
3456       call readi(controlcard,'GENCONSTR',genconstr,0)
3457       call reada(controlcard,'BOXX',boxxsize,100.0d0)
3458       call reada(controlcard,'BOXY',boxysize,100.0d0)
3459       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3460       call readi(controlcard,'TUBEMOD',tubemode,0)
3461       write (iout,*) TUBEmode,"TUBEMODE"
3462       if (TUBEmode.gt.0) then
3463        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3464        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3465        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3466        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3467        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3468        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3469        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3470        buftubebot=bordtubebot+tubebufthick
3471        buftubetop=bordtubetop-tubebufthick
3472       endif
3473
3474 ! CUTOFFF ON ELECTROSTATICS
3475       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3476       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3477       write(iout,*) "R_CUT_ELE=",r_cut_ele
3478 ! Lipidic parameters
3479       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3480       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3481       if (lipthick.gt.0.0d0) then
3482        bordliptop=(boxzsize+lipthick)/2.0
3483        bordlipbot=bordliptop-lipthick
3484       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3485       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3486       buflipbot=bordlipbot+lipbufthick
3487       bufliptop=bordliptop-lipbufthick
3488       if ((lipbufthick*2.0d0).gt.lipthick) &
3489        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3490       endif !lipthick.gt.0
3491       write(iout,*) "bordliptop=",bordliptop
3492       write(iout,*) "bordlipbot=",bordlipbot
3493       write(iout,*) "bufliptop=",bufliptop
3494       write(iout,*) "buflipbot=",buflipbot
3495       write (iout,*) "SHIELD MODE",shield_mode
3496
3497 !C-------------------------
3498       minim=(index(controlcard,'MINIMIZE').gt.0)
3499       dccart=(index(controlcard,'CART').gt.0)
3500       overlapsc=(index(controlcard,'OVERLAP').gt.0)
3501       overlapsc=.not.overlapsc
3502       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3503       searchsc=.not.searchsc
3504       sideadd=(index(controlcard,'SIDEADD').gt.0)
3505       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3506       outpdb=(index(controlcard,'PDBOUT').gt.0)
3507       outmol2=(index(controlcard,'MOL2OUT').gt.0)
3508       pdbref=(index(controlcard,'PDBREF').gt.0)
3509       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3510       indpdb=index(controlcard,'PDBSTART')
3511       extconf=(index(controlcard,'EXTCONF').gt.0)
3512       call readi(controlcard,'IPRINT',iprint,0)
3513       call readi(controlcard,'MAXGEN',maxgen,10000)
3514       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3515       call readi(controlcard,"KDIAG",kdiag,0)
3516       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3517       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3518        write (iout,*) "RESCALE_MODE",rescale_mode
3519       split_ene=index(controlcard,'SPLIT_ENE').gt.0
3520       if (index(controlcard,'REGULAR').gt.0.0D0) then
3521         call reada(controlcard,'WEIDIS',weidis,0.1D0)
3522         modecalc=1
3523         refstr=.true.
3524       endif
3525       if (index(controlcard,'CHECKGRAD').gt.0) then
3526         modecalc=5
3527         if (index(controlcard,'CART').gt.0) then
3528           icheckgrad=1
3529         elseif (index(controlcard,'CARINT').gt.0) then
3530           icheckgrad=2
3531         else
3532           icheckgrad=3
3533         endif
3534       elseif (index(controlcard,'THREAD').gt.0) then
3535         modecalc=2
3536         call readi(controlcard,'THREAD',nthread,0)
3537         if (nthread.gt.0) then
3538           call reada(controlcard,'WEIDIS',weidis,0.1D0)
3539         else
3540           if (fg_rank.eq.0) &
3541           write (iout,'(a)')'A number has to follow the THREAD keyword.'
3542           stop 'Error termination in Read_Control.'
3543         endif
3544       else if (index(controlcard,'MCMA').gt.0) then
3545         modecalc=3
3546       else if (index(controlcard,'MCEE').gt.0) then
3547         modecalc=6
3548       else if (index(controlcard,'MULTCONF').gt.0) then
3549         modecalc=4
3550       else if (index(controlcard,'MAP').gt.0) then
3551         modecalc=7
3552         call readi(controlcard,'MAP',nmap,0)
3553       else if (index(controlcard,'CSA').gt.0) then
3554         modecalc=8
3555 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3556 !rc
3557 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3558 !rc
3559 !rc        modecalc=9
3560 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3561 !fmc        modecalc=10
3562       else if (index(controlcard,'SOFTREG').gt.0) then
3563         modecalc=11
3564       else if (index(controlcard,'CHECK_BOND').gt.0) then
3565         modecalc=-1
3566       else if (index(controlcard,'TEST').gt.0) then
3567         modecalc=-2
3568       else if (index(controlcard,'MD').gt.0) then
3569         modecalc=12
3570       else if (index(controlcard,'RE ').gt.0) then
3571         modecalc=14
3572       endif
3573
3574       lmuca=index(controlcard,'MUCA').gt.0
3575       call readi(controlcard,'MUCADYN',mucadyn,0)
3576       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3577       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3578        then
3579        write (iout,*) 'MUCADYN=',mucadyn
3580        write (iout,*) 'MUCASMOOTH=',muca_smooth
3581       endif
3582
3583       iscode=index(controlcard,'ONE_LETTER')
3584       indphi=index(controlcard,'PHI')
3585       indback=index(controlcard,'BACK')
3586       iranconf=index(controlcard,'RAND_CONF')
3587       i2ndstr=index(controlcard,'USE_SEC_PRED')
3588       gradout=index(controlcard,'GRADOUT').gt.0
3589       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3590       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3591       if (me.eq.king .or. .not.out1file ) &
3592         write (iout,*) "DISTCHAINMAX",distchainmax
3593
3594       if(me.eq.king.or..not.out1file) &
3595        write (iout,'(2a)') diagmeth(kdiag),&
3596         ' routine used to diagonalize matrices.'
3597       if (shield_mode.gt.0) then
3598       pi=3.141592d0
3599 !C VSolvSphere the volume of solving sphere
3600 !C      print *,pi,"pi"
3601 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3602 !C there will be no distinction between proline peptide group and normal peptide
3603 !C group in case of shielding parameters
3604       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3605       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3606       write (iout,*) VSolvSphere,VSolvSphere_div
3607 !C long axis of side chain
3608       do i=1,ntyp
3609       long_r_sidechain(i)=vbldsc0(1,i)
3610       short_r_sidechain(i)=sigma0(i)
3611       write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3612          sigma0(i)
3613       enddo
3614       buff_shield=1.0d0
3615       endif
3616       return
3617       end subroutine read_control
3618 !-----------------------------------------------------------------------------
3619       subroutine read_REMDpar
3620 !
3621 ! Read REMD settings
3622 !
3623 !       use control
3624 !       use energy
3625 !       use geometry
3626       use REMD_data
3627       use MPI_data
3628       use control_data, only:out1file
3629 !      implicit real*8 (a-h,o-z)
3630 !      include 'DIMENSIONS'
3631 !      include 'COMMON.IOUNITS'
3632 !      include 'COMMON.TIME1'
3633 !      include 'COMMON.MD'
3634       use MD_data
3635 !el #ifndef LANG0
3636 !el      include 'COMMON.LANGEVIN'
3637 !el #else
3638 !el      include 'COMMON.LANGEVIN.lang0'
3639 !el #endif
3640 !      include 'COMMON.INTERACT'
3641 !      include 'COMMON.NAMES'
3642 !      include 'COMMON.GEO'
3643 !      include 'COMMON.REMD'
3644 !      include 'COMMON.CONTROL'
3645 !      include 'COMMON.SETUP'
3646 !      character(len=80) :: ucase
3647       character(len=320) :: controlcard
3648       character(len=3200) :: controlcard1
3649       integer :: iremd_m_total
3650 !el local variables
3651       integer :: i
3652 !     real(kind=8) :: var,ene
3653
3654       if(me.eq.king.or..not.out1file) &
3655        write (iout,*) "REMD setup"
3656
3657       call card_concat(controlcard,.true.)
3658       call readi(controlcard,"NREP",nrep,3)
3659       call readi(controlcard,"NSTEX",nstex,1000)
3660       call reada(controlcard,"RETMIN",retmin,10.0d0)
3661       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3662       mremdsync=(index(controlcard,'SYNC').gt.0)
3663       call readi(controlcard,"NSYN",i_sync_step,100)
3664       restart1file=(index(controlcard,'REST1FILE').gt.0)
3665       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3666       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3667       if(max_cache_traj_use.gt.max_cache_traj) &
3668                 max_cache_traj_use=max_cache_traj
3669       if(me.eq.king.or..not.out1file) then
3670 !d       if (traj1file) then
3671 !rc caching is in testing - NTWX is not ignored
3672 !d        write (iout,*) "NTWX value is ignored"
3673 !d        write (iout,*) "  trajectory is stored to one file by master"
3674 !d        write (iout,*) "  before exchange at NSTEX intervals"
3675 !d       endif
3676        write (iout,*) "NREP= ",nrep
3677        write (iout,*) "NSTEX= ",nstex
3678        write (iout,*) "SYNC= ",mremdsync
3679        write (iout,*) "NSYN= ",i_sync_step
3680        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3681       endif
3682       remd_tlist=.false.
3683       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3684       if (index(controlcard,'TLIST').gt.0) then
3685          remd_tlist=.true.
3686          call card_concat(controlcard1,.true.)
3687          read(controlcard1,*) (remd_t(i),i=1,nrep)
3688          if(me.eq.king.or..not.out1file) &
3689           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3690       endif
3691       remd_mlist=.false.
3692       if (index(controlcard,'MLIST').gt.0) then
3693          remd_mlist=.true.
3694          call card_concat(controlcard1,.true.)
3695          read(controlcard1,*) (remd_m(i),i=1,nrep)
3696          if(me.eq.king.or..not.out1file) then
3697           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3698           iremd_m_total=0
3699           do i=1,nrep
3700            iremd_m_total=iremd_m_total+remd_m(i)
3701           enddo
3702           write (iout,*) 'Total number of replicas ',iremd_m_total
3703          endif
3704       endif
3705       if(me.eq.king.or..not.out1file) &
3706        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3707       return
3708       end subroutine read_REMDpar
3709 !-----------------------------------------------------------------------------
3710       subroutine read_MDpar
3711 !
3712 ! Read MD settings
3713 !
3714       use control_data, only: r_cut,rlamb,out1file
3715       use energy_data
3716       use geometry_data, only: pi
3717       use MPI_data
3718 !      implicit real*8 (a-h,o-z)
3719 !      include 'DIMENSIONS'
3720 !      include 'COMMON.IOUNITS'
3721 !      include 'COMMON.TIME1'
3722 !      include 'COMMON.MD'
3723       use MD_data
3724 !el #ifndef LANG0
3725 !el      include 'COMMON.LANGEVIN'
3726 !el #else
3727 !el      include 'COMMON.LANGEVIN.lang0'
3728 !el #endif
3729 !      include 'COMMON.INTERACT'
3730 !      include 'COMMON.NAMES'
3731 !      include 'COMMON.GEO'
3732 !      include 'COMMON.SETUP'
3733 !      include 'COMMON.CONTROL'
3734 !      include 'COMMON.SPLITELE'
3735 !      character(len=80) :: ucase
3736       character(len=320) :: controlcard
3737 !el local variables
3738       integer :: i,j
3739       real(kind=8) :: eta
3740
3741       call card_concat(controlcard,.true.)
3742       call readi(controlcard,"NSTEP",n_timestep,1000000)
3743       call readi(controlcard,"NTWE",ntwe,100)
3744       call readi(controlcard,"NTWX",ntwx,1000)
3745       call reada(controlcard,"DT",d_time,1.0d-1)
3746       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3747       call reada(controlcard,"DAMAX",damax,1.0d1)
3748       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3749       call readi(controlcard,"LANG",lang,0)
3750       RESPA = index(controlcard,"RESPA") .gt. 0
3751       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3752       ntime_split0=ntime_split
3753       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3754       ntime_split0=ntime_split
3755       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3756       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3757       rest = index(controlcard,"REST").gt.0
3758       tbf = index(controlcard,"TBF").gt.0
3759       usampl = index(controlcard,"USAMPL").gt.0
3760       mdpdb = index(controlcard,"MDPDB").gt.0
3761       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3762       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3763       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3764       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3765       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3766       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3767       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3768       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3769       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3770       large = index(controlcard,"LARGE").gt.0
3771       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3772       rattle = index(controlcard,"RATTLE").gt.0
3773 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3774       nset=0
3775       if(usampl) then
3776         call read_fragments
3777       endif
3778
3779       if(me.eq.king.or..not.out1file) then
3780        write (iout,*)
3781        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3782        write (iout,*)
3783        write (iout,'(a)') "The units are:"
3784        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3785        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3786         " acceleration: angstrom/(48.9 fs)**2"
3787        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3788        write (iout,*)
3789        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3790        write (iout,'(a60,f10.5,a)') &
3791         "Initial time step of numerical integration:",d_time,&
3792         " natural units"
3793        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3794        if (RESPA) then
3795         write (iout,'(2a,i4,a)') &
3796           "A-MTS algorithm used; initial time step for fast-varying",&
3797           " short-range forces split into",ntime_split," steps."
3798         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3799          r_cut," lambda",rlamb
3800        endif
3801        write (iout,'(2a,f10.5)') &
3802         "Maximum acceleration threshold to reduce the time step",&
3803         "/increase split number:",damax
3804        write (iout,'(2a,f10.5)') &
3805         "Maximum predicted energy drift to reduce the timestep",&
3806         "/increase split number:",edriftmax
3807        write (iout,'(a60,f10.5)') &
3808        "Maximum velocity threshold to reduce velocities:",dvmax
3809        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3810        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3811        if (rattle) write (iout,'(a60)') &
3812         "Rattle algorithm used to constrain the virtual bonds"
3813       endif
3814       reset_fricmat=1000
3815       if (lang.gt.0) then
3816         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3817         call reada(controlcard,"RWAT",rwat,1.4d0)
3818         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3819         surfarea=index(controlcard,"SURFAREA").gt.0
3820         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3821         if(me.eq.king.or..not.out1file)then
3822          write (iout,'(/a,$)') "Langevin dynamics calculation"
3823          if (lang.eq.1) then
3824           write (iout,'(a/)') &
3825             " with direct integration of Langevin equations"
3826          else if (lang.eq.2) then
3827           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3828          else if (lang.eq.3) then
3829           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3830          else if (lang.eq.4) then
3831           write (iout,'(a/)') " in overdamped mode"
3832          else
3833           write (iout,'(//a,i5)') &
3834             "=========== ERROR: Unknown Langevin dynamics mode:",lang
3835           stop
3836          endif
3837          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3838          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3839          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3840          write (iout,'(a60,f10.5)') &
3841          "Scaling factor of the friction forces:",scal_fric
3842          if (surfarea) write (iout,'(2a,i10,a)') &
3843            "Friction coefficients will be scaled by solvent-accessible",&
3844            " surface area every",reset_fricmat," steps."
3845         endif
3846 ! Calculate friction coefficients and bounds of stochastic forces
3847         eta=6*pi*cPoise*etawat
3848         if(me.eq.king.or..not.out1file) &
3849          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3850           eta
3851 !        allocate(gamp
3852         do j=1,5 !types of molecules
3853         gamp(j)=scal_fric*(pstok(j)+rwat)*eta
3854         stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
3855         enddo
3856         allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
3857         do j=1,5 !types of molecules
3858         do i=1,ntyp
3859           gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
3860           stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
3861         enddo
3862         enddo
3863
3864         if(me.eq.king.or..not.out1file)then
3865          write (iout,'(/2a/)') &
3866          "Radii of site types and friction coefficients and std's of",&
3867          " stochastic forces of fully exposed sites"
3868          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
3869          do i=1,ntyp
3870           write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
3871            gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
3872          enddo
3873         endif
3874       else if (tbf) then
3875         if(me.eq.king.or..not.out1file)then
3876          write (iout,'(a)') "Berendsen bath calculation"
3877          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3878          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3879          if (reset_moment) &
3880          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3881          count_reset_moment," steps"
3882          if (reset_vel) &
3883           write (iout,'(a,i10,a)') &
3884           "Velocities will be reset at random every",count_reset_vel,&
3885          " steps"
3886         endif
3887       else
3888         if(me.eq.king.or..not.out1file) &
3889          write (iout,'(a31)') "Microcanonical mode calculation"
3890       endif
3891       if(me.eq.king.or..not.out1file)then
3892        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3893        if (usampl) then
3894           write(iout,*) "MD running with constraints."
3895           write(iout,*) "Equilibration time ", eq_time, " mtus."
3896           write(iout,*) "Constraining ", nfrag," fragments."
3897           write(iout,*) "Length of each fragment, weight and q0:"
3898           do iset=1,nset
3899            write (iout,*) "Set of restraints #",iset
3900            do i=1,nfrag
3901               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3902                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3903            enddo
3904            write(iout,*) "constraints between ", npair, "fragments."
3905            write(iout,*) "constraint pairs, weights and q0:"
3906            do i=1,npair
3907             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3908                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3909            enddo
3910            write(iout,*) "angle constraints within ", nfrag_back,&
3911             "backbone fragments."
3912            write(iout,*) "fragment, weights:"
3913            do i=1,nfrag_back
3914             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3915                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3916                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3917            enddo
3918           enddo
3919         iset=mod(kolor,nset)+1
3920        endif
3921       endif
3922       if(me.eq.king.or..not.out1file) &
3923        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3924       return
3925       end subroutine read_MDpar
3926 !-----------------------------------------------------------------------------
3927       subroutine map_read
3928
3929       use map_data
3930 !      implicit real*8 (a-h,o-z)
3931 !      include 'DIMENSIONS'
3932 !      include 'COMMON.MAP'
3933 !      include 'COMMON.IOUNITS'
3934       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3935       character(len=80) :: mapcard      !,ucase
3936 !el local variables
3937       integer :: imap
3938 !     real(kind=8) :: var,ene
3939
3940       do imap=1,nmap
3941         read (inp,'(a)') mapcard
3942         mapcard=ucase(mapcard)
3943         if (index(mapcard,'PHI').gt.0) then
3944           kang(imap)=1
3945         else if (index(mapcard,'THE').gt.0) then
3946           kang(imap)=2
3947         else if (index(mapcard,'ALP').gt.0) then
3948           kang(imap)=3
3949         else if (index(mapcard,'OME').gt.0) then
3950           kang(imap)=4
3951         else
3952           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3953           stop 'Error - illegal variable spec in MAP card.'
3954         endif
3955         call readi (mapcard,'RES1',res1(imap),0)
3956         call readi (mapcard,'RES2',res2(imap),0)
3957         if (res1(imap).eq.0) then
3958           res1(imap)=res2(imap)
3959         else if (res2(imap).eq.0) then
3960           res2(imap)=res1(imap)
3961         endif
3962         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3963           write (iout,'(a)') &
3964           'Error - illegal definition of variable group in MAP.'
3965           stop 'Error - illegal definition of variable group in MAP.'
3966         endif
3967         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3968         call reada(mapcard,'TO',ang_to(imap),0.0D0)
3969         call readi(mapcard,'NSTEP',nstep(imap),0)
3970         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3971           write (iout,'(a)') &
3972            'Illegal boundary and/or step size specification in MAP.'
3973           stop 'Illegal boundary and/or step size specification in MAP.'
3974         endif
3975       enddo ! imap
3976       return
3977       end subroutine map_read
3978 !-----------------------------------------------------------------------------
3979       subroutine csaread
3980
3981       use control_data, only: vdisulf
3982       use csa_data
3983 !      implicit real*8 (a-h,o-z)
3984 !      include 'DIMENSIONS'
3985 !      include 'COMMON.IOUNITS'
3986 !      include 'COMMON.GEO'
3987 !      include 'COMMON.CSA'
3988 !      include 'COMMON.BANK'
3989 !      include 'COMMON.CONTROL'
3990 !      character(len=80) :: ucase
3991       character(len=620) :: mcmcard
3992 !el local variables
3993 !     integer :: ntf,ik,iw_pdb
3994 !     real(kind=8) :: var,ene
3995
3996       call card_concat(mcmcard,.true.)
3997
3998       call readi(mcmcard,'NCONF',nconf,50)
3999       call readi(mcmcard,'NADD',nadd,0)
4000       call readi(mcmcard,'JSTART',jstart,1)
4001       call readi(mcmcard,'JEND',jend,1)
4002       call readi(mcmcard,'NSTMAX',nstmax,500000)
4003       call readi(mcmcard,'N0',n0,1)
4004       call readi(mcmcard,'N1',n1,6)
4005       call readi(mcmcard,'N2',n2,4)
4006       call readi(mcmcard,'N3',n3,0)
4007       call readi(mcmcard,'N4',n4,0)
4008       call readi(mcmcard,'N5',n5,0)
4009       call readi(mcmcard,'N6',n6,10)
4010       call readi(mcmcard,'N7',n7,0)
4011       call readi(mcmcard,'N8',n8,0)
4012       call readi(mcmcard,'N9',n9,0)
4013       call readi(mcmcard,'N14',n14,0)
4014       call readi(mcmcard,'N15',n15,0)
4015       call readi(mcmcard,'N16',n16,0)
4016       call readi(mcmcard,'N17',n17,0)
4017       call readi(mcmcard,'N18',n18,0)
4018
4019       vdisulf=(index(mcmcard,'DYNSS').gt.0)
4020
4021       call readi(mcmcard,'NDIFF',ndiff,2)
4022       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4023       call readi(mcmcard,'IS1',is1,1)
4024       call readi(mcmcard,'IS2',is2,8)
4025       call readi(mcmcard,'NRAN0',nran0,4)
4026       call readi(mcmcard,'NRAN1',nran1,2)
4027       call readi(mcmcard,'IRR',irr,1)
4028       call readi(mcmcard,'NSEED',nseed,20)
4029       call readi(mcmcard,'NTOTAL',ntotal,10000)
4030       call reada(mcmcard,'CUT1',cut1,2.0d0)
4031       call reada(mcmcard,'CUT2',cut2,5.0d0)
4032       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4033       call readi(mcmcard,'ICMAX',icmax,3)
4034       call readi(mcmcard,'IRESTART',irestart,0)
4035 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
4036       ntbankm=0
4037 !!bankt
4038       call reada(mcmcard,'DELE',dele,20.0d0)
4039       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4040       call readi(mcmcard,'IREF',iref,0)
4041       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4042       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4043       call readi(mcmcard,'NCONF_IN',nconf_in,0)
4044       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4045       write (iout,*) "NCONF_IN",nconf_in
4046       return
4047       end subroutine csaread
4048 !-----------------------------------------------------------------------------
4049       subroutine mcmread
4050
4051       use mcm_data
4052       use control_data, only: MaxMoveType
4053       use MD_data
4054       use minim_data
4055 !      implicit real*8 (a-h,o-z)
4056 !      include 'DIMENSIONS'
4057 !      include 'COMMON.MCM'
4058 !      include 'COMMON.MCE'
4059 !      include 'COMMON.IOUNITS'
4060 !      character(len=80) :: ucase
4061       character(len=320) :: mcmcard
4062 !el local variables
4063       integer :: i
4064 !     real(kind=8) :: var,ene
4065
4066       call card_concat(mcmcard,.true.)
4067       call readi(mcmcard,'MAXACC',maxacc,100)
4068       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4069       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4070       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4071       call readi(mcmcard,'MAXREPM',maxrepm,200)
4072       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4073       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4074       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4075       call reada(mcmcard,'E_UP',e_up,5.0D0)
4076       call reada(mcmcard,'DELTE',delte,0.1D0)
4077       call readi(mcmcard,'NSWEEP',nsweep,5)
4078       call readi(mcmcard,'NSTEPH',nsteph,0)
4079       call readi(mcmcard,'NSTEPC',nstepc,0)
4080       call reada(mcmcard,'TMIN',tmin,298.0D0)
4081       call reada(mcmcard,'TMAX',tmax,298.0D0)
4082       call readi(mcmcard,'NWINDOW',nwindow,0)
4083       call readi(mcmcard,'PRINT_MC',print_mc,0)
4084       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4085       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4086       ent_read=(index(mcmcard,'ENT_READ').gt.0)
4087       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4088       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4089       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4090       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4091       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4092       if (nwindow.gt.0) then
4093         allocate(winstart(nwindow))     !!el (maxres)
4094         allocate(winend(nwindow))       !!el
4095         allocate(winlen(nwindow))       !!el
4096         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4097         do i=1,nwindow
4098           winlen(i)=winend(i)-winstart(i)+1
4099         enddo
4100       endif
4101       if (tmax.lt.tmin) tmax=tmin
4102       if (tmax.eq.tmin) then
4103         nstepc=0
4104         nsteph=0
4105       endif
4106       if (nstepc.gt.0 .and. nsteph.gt.0) then
4107         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4108         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4109       endif
4110       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4111 ! Probabilities of different move types
4112       sumpro_type(0)=0.0D0
4113       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4114       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4115       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4116       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
4117       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4118       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4119       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4120       do i=1,MaxMoveType
4121         print *,'i',i,' sumprotype',sumpro_type(i)
4122         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4123         print *,'i',i,' sumprotype',sumpro_type(i)
4124       enddo
4125       return
4126       end subroutine mcmread
4127 !-----------------------------------------------------------------------------
4128       subroutine read_minim
4129
4130       use minim_data
4131 !      implicit real*8 (a-h,o-z)
4132 !      include 'DIMENSIONS'
4133 !      include 'COMMON.MINIM'
4134 !      include 'COMMON.IOUNITS'
4135 !      character(len=80) :: ucase
4136       character(len=320) :: minimcard
4137 !el local variables
4138 !     integer :: ntf,ik,iw_pdb
4139 !     real(kind=8) :: var,ene
4140
4141       call card_concat(minimcard,.true.)
4142       call readi(minimcard,'MAXMIN',maxmin,2000)
4143       call readi(minimcard,'MAXFUN',maxfun,5000)
4144       call readi(minimcard,'MINMIN',minmin,maxmin)
4145       call readi(minimcard,'MINFUN',minfun,maxmin)
4146       call reada(minimcard,'TOLF',tolf,1.0D-2)
4147       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4148       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4149       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4150       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4151       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4152                'Options in energy minimization:'
4153       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4154        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4155        'MinMin:',MinMin,' MinFun:',MinFun,&
4156        ' TolF:',TolF,' RTolF:',RTolF
4157       return
4158       end subroutine read_minim
4159 !-----------------------------------------------------------------------------
4160       subroutine openunits
4161
4162       use MD_data, only: usampl
4163       use csa_data
4164       use MPI_data
4165       use control_data, only:out1file
4166       use control, only: getenv_loc
4167 !      implicit real*8 (a-h,o-z)
4168 !      include 'DIMENSIONS'
4169 #ifdef MPI
4170       include 'mpif.h'
4171       character(len=16) :: form,nodename
4172       integer :: nodelen,ierror,npos
4173 #endif
4174 !      include 'COMMON.SETUP'
4175 !      include 'COMMON.IOUNITS'
4176 !      include 'COMMON.MD'
4177 !      include 'COMMON.CONTROL'
4178       integer :: lenpre,lenpot,lentmp   !,ilen
4179 !el      external ilen
4180       character(len=3) :: out1file_text !,ucase
4181       character(len=3) :: ll
4182 !el      external ucase
4183 !el local variables
4184 !     integer :: ntf,ik,iw_pdb
4185 !     real(kind=8) :: var,ene
4186 !
4187 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4188       call getenv_loc("PREFIX",prefix)
4189       pref_orig = prefix
4190       call getenv_loc("POT",pot)
4191       call getenv_loc("DIRTMP",tmpdir)
4192       call getenv_loc("CURDIR",curdir)
4193       call getenv_loc("OUT1FILE",out1file_text)
4194 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4195       out1file_text=ucase(out1file_text)
4196       if (out1file_text(1:1).eq."Y") then
4197         out1file=.true.
4198       else
4199         out1file=fg_rank.gt.0
4200       endif
4201       lenpre=ilen(prefix)
4202       lenpot=ilen(pot)
4203       lentmp=ilen(tmpdir)
4204       if (lentmp.gt.0) then
4205           write (*,'(80(1h!))')
4206           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
4207           write (*,'(80(1h!))')
4208           write (*,*)"All output files will be on node /tmp directory."
4209 #ifdef MPI
4210         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4211         if (me.eq.king) then
4212           write (*,*) "The master node is ",nodename
4213         else if (fg_rank.eq.0) then
4214           write (*,*) "I am the CG slave node ",nodename
4215         else
4216           write (*,*) "I am the FG slave node ",nodename
4217         endif
4218 #endif
4219         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4220         lenpre = lentmp+lenpre+1
4221       endif
4222       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4223 ! Get the names and open the input files
4224 #if defined(WINIFL) || defined(WINPGI)
4225       open(1,file=pref_orig(:ilen(pref_orig))// &
4226         '.inp',status='old',readonly,shared)
4227        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4228 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4229 ! Get parameter filenames and open the parameter files.
4230       call getenv_loc('BONDPAR',bondname)
4231       open (ibond,file=bondname,status='old',readonly,shared)
4232       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4233       open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4234       call getenv_loc('THETPAR',thetname)
4235       open (ithep,file=thetname,status='old',readonly,shared)
4236       call getenv_loc('ROTPAR',rotname)
4237       open (irotam,file=rotname,status='old',readonly,shared)
4238       call getenv_loc('TORPAR',torname)
4239       open (itorp,file=torname,status='old',readonly,shared)
4240       call getenv_loc('TORDPAR',tordname)
4241       open (itordp,file=tordname,status='old',readonly,shared)
4242       call getenv_loc('FOURIER',fouriername)
4243       open (ifourier,file=fouriername,status='old',readonly,shared)
4244       call getenv_loc('ELEPAR',elename)
4245       open (ielep,file=elename,status='old',readonly,shared)
4246       call getenv_loc('SIDEPAR',sidename)
4247       open (isidep,file=sidename,status='old',readonly,shared)
4248
4249       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4250       open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4251       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4252       open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4253       call getenv_loc('TORPAR_NUCL',torname_nucl)
4254       open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4255       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4256       open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4257       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4258       open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4259
4260
4261 #elif (defined CRAY) || (defined AIX)
4262       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4263         action='read')
4264 !      print *,"Processor",myrank," opened file 1"
4265       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4266 !      print *,"Processor",myrank," opened file 9"
4267 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4268 ! Get parameter filenames and open the parameter files.
4269       call getenv_loc('BONDPAR',bondname)
4270       open (ibond,file=bondname,status='old',action='read')
4271       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4272       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4273
4274 !      print *,"Processor",myrank," opened file IBOND"
4275       call getenv_loc('THETPAR',thetname)
4276       open (ithep,file=thetname,status='old',action='read')
4277 !      print *,"Processor",myrank," opened file ITHEP"
4278       call getenv_loc('ROTPAR',rotname)
4279       open (irotam,file=rotname,status='old',action='read')
4280 !      print *,"Processor",myrank," opened file IROTAM"
4281       call getenv_loc('TORPAR',torname)
4282       open (itorp,file=torname,status='old',action='read')
4283 !      print *,"Processor",myrank," opened file ITORP"
4284       call getenv_loc('TORDPAR',tordname)
4285       open (itordp,file=tordname,status='old',action='read')
4286 !      print *,"Processor",myrank," opened file ITORDP"
4287       call getenv_loc('SCCORPAR',sccorname)
4288       open (isccor,file=sccorname,status='old',action='read')
4289 !      print *,"Processor",myrank," opened file ISCCOR"
4290       call getenv_loc('FOURIER',fouriername)
4291       open (ifourier,file=fouriername,status='old',action='read')
4292 !      print *,"Processor",myrank," opened file IFOURIER"
4293       call getenv_loc('ELEPAR',elename)
4294       open (ielep,file=elename,status='old',action='read')
4295 !      print *,"Processor",myrank," opened file IELEP"
4296       call getenv_loc('SIDEPAR',sidename)
4297       open (isidep,file=sidename,status='old',action='read')
4298
4299       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4300       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4301       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4302       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4303       call getenv_loc('TORPAR_NUCL',torname_nucl)
4304       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4305       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4306       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4307       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4308       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4309
4310       call getenv_loc('LIPTRANPAR',liptranname)
4311       open (iliptranpar,file=liptranname,status='old',action='read')
4312       call getenv_loc('TUBEPAR',tubename)
4313       open (itube,file=tubename,status='old',action='read')
4314       call getenv_loc('IONPAR',ionname)
4315       open (iion,file=ionname,status='old',action='read')
4316
4317 !      print *,"Processor",myrank," opened file ISIDEP"
4318 !      print *,"Processor",myrank," opened parameter files"
4319 #elif (defined G77)
4320       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4321       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4322 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4323 ! Get parameter filenames and open the parameter files.
4324       call getenv_loc('BONDPAR',bondname)
4325       open (ibond,file=bondname,status='old')
4326       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4327       open (ibond_nucl,file=bondname_nucl,status='old')
4328
4329       call getenv_loc('THETPAR',thetname)
4330       open (ithep,file=thetname,status='old')
4331       call getenv_loc('ROTPAR',rotname)
4332       open (irotam,file=rotname,status='old')
4333       call getenv_loc('TORPAR',torname)
4334       open (itorp,file=torname,status='old')
4335       call getenv_loc('TORDPAR',tordname)
4336       open (itordp,file=tordname,status='old')
4337       call getenv_loc('SCCORPAR',sccorname)
4338       open (isccor,file=sccorname,status='old')
4339       call getenv_loc('FOURIER',fouriername)
4340       open (ifourier,file=fouriername,status='old')
4341       call getenv_loc('ELEPAR',elename)
4342       open (ielep,file=elename,status='old')
4343       call getenv_loc('SIDEPAR',sidename)
4344       open (isidep,file=sidename,status='old')
4345
4346       open (ithep_nucl,file=thetname_nucl,status='old')
4347       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4348       open (irotam_nucl,file=rotname_nucl,status='old')
4349       call getenv_loc('TORPAR_NUCL',torname_nucl)
4350       open (itorp_nucl,file=torname_nucl,status='old')
4351       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4352       open (itordp_nucl,file=tordname_nucl,status='old')
4353       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4354       open (isidep_nucl,file=sidename_nucl,status='old')
4355
4356       call getenv_loc('LIPTRANPAR',liptranname)
4357       open (iliptranpar,file=liptranname,status='old')
4358       call getenv_loc('TUBEPAR',tubename)
4359       open (itube,file=tubename,status='old')
4360       call getenv_loc('IONPAR',ionname)
4361       open (iion,file=ionname,status='old')
4362 #else
4363       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4364         readonly)
4365        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4366 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4367 ! Get parameter filenames and open the parameter files.
4368       call getenv_loc('BONDPAR',bondname)
4369       open (ibond,file=bondname,status='old',action='read')
4370       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4371       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4372       call getenv_loc('THETPAR',thetname)
4373       open (ithep,file=thetname,status='old',action='read')
4374       call getenv_loc('ROTPAR',rotname)
4375       open (irotam,file=rotname,status='old',action='read')
4376       call getenv_loc('TORPAR',torname)
4377       open (itorp,file=torname,status='old',action='read')
4378       call getenv_loc('TORDPAR',tordname)
4379       open (itordp,file=tordname,status='old',action='read')
4380       call getenv_loc('SCCORPAR',sccorname)
4381       open (isccor,file=sccorname,status='old',action='read')
4382 #ifndef CRYST_THETA
4383       call getenv_loc('THETPARPDB',thetname_pdb)
4384       print *,"thetname_pdb ",thetname_pdb
4385       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4386       print *,ithep_pdb," opened"
4387 #endif
4388       call getenv_loc('FOURIER',fouriername)
4389       open (ifourier,file=fouriername,status='old',readonly)
4390       call getenv_loc('ELEPAR',elename)
4391       open (ielep,file=elename,status='old',readonly)
4392       call getenv_loc('SIDEPAR',sidename)
4393       open (isidep,file=sidename,status='old',readonly)
4394
4395       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4396       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4397       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4398       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4399       call getenv_loc('TORPAR_NUCL',torname_nucl)
4400       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4401       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4402       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4403       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4404       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4405
4406       call getenv_loc('LIPTRANPAR',liptranname)
4407       open (iliptranpar,file=liptranname,status='old',action='read')
4408       call getenv_loc('TUBEPAR',tubename)
4409       open (itube,file=tubename,status='old',action='read')
4410       call getenv_loc('IONPAR',ionname)
4411       open (iion,file=ionname,status='old',action='read')
4412
4413 #ifndef CRYST_SC
4414       call getenv_loc('ROTPARPDB',rotname_pdb)
4415       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4416 #endif
4417 #endif
4418       call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4419 #if defined(WINIFL) || defined(WINPGI)
4420       open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4421 #elif (defined CRAY)  || (defined AIX)
4422       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4423 #elif (defined G77)
4424       open (iscpp_nucl,file=scpname_nucl,status='old')
4425 #else
4426       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4427 #endif
4428
4429 #ifndef OLDSCP
4430 !
4431 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4432 ! Use -DOLDSCP to use hard-coded constants instead.
4433 !
4434       call getenv_loc('SCPPAR',scpname)
4435 #if defined(WINIFL) || defined(WINPGI)
4436       open (iscpp,file=scpname,status='old',readonly,shared)
4437 #elif (defined CRAY)  || (defined AIX)
4438       open (iscpp,file=scpname,status='old',action='read')
4439 #elif (defined G77)
4440       open (iscpp,file=scpname,status='old')
4441 #else
4442       open (iscpp,file=scpname,status='old',action='read')
4443 #endif
4444 #endif
4445       call getenv_loc('PATTERN',patname)
4446 #if defined(WINIFL) || defined(WINPGI)
4447       open (icbase,file=patname,status='old',readonly,shared)
4448 #elif (defined CRAY)  || (defined AIX)
4449       open (icbase,file=patname,status='old',action='read')
4450 #elif (defined G77)
4451       open (icbase,file=patname,status='old')
4452 #else
4453       open (icbase,file=patname,status='old',action='read')
4454 #endif
4455 #ifdef MPI
4456 ! Open output file only for CG processes
4457 !      print *,"Processor",myrank," fg_rank",fg_rank
4458       if (fg_rank.eq.0) then
4459
4460       if (nodes.eq.1) then
4461         npos=3
4462       else
4463         npos = dlog10(dfloat(nodes-1))+1
4464       endif
4465       if (npos.lt.3) npos=3
4466       write (liczba,'(i1)') npos
4467       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4468         //')'
4469       write (liczba,form) me
4470       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4471         liczba(:ilen(liczba))
4472       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4473         //'.int'
4474       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4475         //'.pdb'
4476       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4477         liczba(:ilen(liczba))//'.mol2'
4478       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4479         liczba(:ilen(liczba))//'.stat'
4480       if (lentmp.gt.0) &
4481         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4482             //liczba(:ilen(liczba))//'.stat')
4483       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4484         //'.rst'
4485       if(usampl) then
4486           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4487        liczba(:ilen(liczba))//'.const'
4488       endif
4489
4490       endif
4491 #else
4492       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4493       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4494       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4495       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4496       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4497       if (lentmp.gt.0) &
4498         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4499           '.stat')
4500       rest2name=prefix(:ilen(prefix))//'.rst'
4501       if(usampl) then
4502          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4503       endif
4504 #endif
4505 #if defined(AIX) || defined(PGI)
4506       if (me.eq.king .or. .not. out1file) &
4507          open(iout,file=outname,status='unknown')
4508 #ifdef DEBUG
4509       if (fg_rank.gt.0) then
4510         write (liczba,'(i3.3)') myrank/nfgtasks
4511         write (ll,'(bz,i3.3)') fg_rank
4512         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4513          status='unknown')
4514       endif
4515 #endif
4516       if(me.eq.king) then
4517        open(igeom,file=intname,status='unknown',position='append')
4518        open(ipdb,file=pdbname,status='unknown')
4519        open(imol2,file=mol2name,status='unknown')
4520        open(istat,file=statname,status='unknown',position='append')
4521       else
4522 !1out       open(iout,file=outname,status='unknown')
4523       endif
4524 #else
4525       if (me.eq.king .or. .not.out1file) &
4526           open(iout,file=outname,status='unknown')
4527 #ifdef DEBUG
4528       if (fg_rank.gt.0) then
4529         write (liczba,'(i3.3)') myrank/nfgtasks
4530         write (ll,'(bz,i3.3)') fg_rank
4531         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4532          status='unknown')
4533       endif
4534 #endif
4535       if(me.eq.king) then
4536        open(igeom,file=intname,status='unknown',access='append')
4537        open(ipdb,file=pdbname,status='unknown')
4538        open(imol2,file=mol2name,status='unknown')
4539        open(istat,file=statname,status='unknown',access='append')
4540       else
4541 !1out       open(iout,file=outname,status='unknown')
4542       endif
4543 #endif
4544       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4545       csa_seed=prefix(:lenpre)//'.CSA.seed'
4546       csa_history=prefix(:lenpre)//'.CSA.history'
4547       csa_bank=prefix(:lenpre)//'.CSA.bank'
4548       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4549       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4550       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4551 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4552       csa_int=prefix(:lenpre)//'.int'
4553       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4554       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4555       csa_in=prefix(:lenpre)//'.CSA.in'
4556 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4557 ! Write file names
4558       if (me.eq.king)then
4559       write (iout,'(80(1h-))')
4560       write (iout,'(30x,a)') "FILE ASSIGNMENT"
4561       write (iout,'(80(1h-))')
4562       write (iout,*) "Input file                      : ",&
4563         pref_orig(:ilen(pref_orig))//'.inp'
4564       write (iout,*) "Output file                     : ",&
4565         outname(:ilen(outname))
4566       write (iout,*)
4567       write (iout,*) "Sidechain potential file        : ",&
4568         sidename(:ilen(sidename))
4569 #ifndef OLDSCP
4570       write (iout,*) "SCp potential file              : ",&
4571         scpname(:ilen(scpname))
4572 #endif
4573       write (iout,*) "Electrostatic potential file    : ",&
4574         elename(:ilen(elename))
4575       write (iout,*) "Cumulant coefficient file       : ",&
4576         fouriername(:ilen(fouriername))
4577       write (iout,*) "Torsional parameter file        : ",&
4578         torname(:ilen(torname))
4579       write (iout,*) "Double torsional parameter file : ",&
4580         tordname(:ilen(tordname))
4581       write (iout,*) "SCCOR parameter file : ",&
4582         sccorname(:ilen(sccorname))
4583       write (iout,*) "Bond & inertia constant file    : ",&
4584         bondname(:ilen(bondname))
4585       write (iout,*) "Bending parameter file          : ",&
4586         thetname(:ilen(thetname))
4587       write (iout,*) "Rotamer parameter file          : ",&
4588         rotname(:ilen(rotname))
4589 !el----
4590 #ifndef CRYST_THETA
4591       write (iout,*) "Thetpdb parameter file          : ",&
4592         thetname_pdb(:ilen(thetname_pdb))
4593 #endif
4594 !el
4595       write (iout,*) "Threading database              : ",&
4596         patname(:ilen(patname))
4597       if (lentmp.ne.0) &
4598       write (iout,*)" DIRTMP                          : ",&
4599         tmpdir(:lentmp)
4600       write (iout,'(80(1h-))')
4601       endif
4602       return
4603       end subroutine openunits
4604 !-----------------------------------------------------------------------------
4605       subroutine readrst
4606
4607       use geometry_data, only: nres,dc
4608       use MD_data
4609 !      implicit real*8 (a-h,o-z)
4610 !      include 'DIMENSIONS'
4611 !      include 'COMMON.CHAIN'
4612 !      include 'COMMON.IOUNITS'
4613 !      include 'COMMON.MD'
4614 !el local variables
4615       integer ::i,j
4616 !     real(kind=8) :: var,ene
4617
4618       open(irest2,file=rest2name,status='unknown')
4619       read(irest2,*) totT,EK,potE,totE,t_bath
4620       totTafm=totT
4621 !      do i=1,2*nres
4622 ! AL 4/17/17: Now reading d_t(0,:) too
4623       do i=0,2*nres
4624          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4625       enddo
4626 !      do i=1,2*nres
4627 ! AL 4/17/17: Now reading d_c(0,:) too
4628       do i=0,2*nres
4629          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4630       enddo
4631       if(usampl) then
4632              read (irest2,*) iset
4633       endif
4634       close(irest2)
4635       return
4636       end subroutine readrst
4637 !-----------------------------------------------------------------------------
4638       subroutine read_fragments
4639
4640       use energy_data
4641 !      use geometry
4642       use control_data, only:out1file
4643       use MD_data
4644       use MPI_data
4645 !      implicit real*8 (a-h,o-z)
4646 !      include 'DIMENSIONS'
4647 #ifdef MPI
4648       include 'mpif.h'
4649 #endif
4650 !      include 'COMMON.SETUP'
4651 !      include 'COMMON.CHAIN'
4652 !      include 'COMMON.IOUNITS'
4653 !      include 'COMMON.MD'
4654 !      include 'COMMON.CONTROL'
4655 !el local variables
4656       integer :: i
4657 !     real(kind=8) :: var,ene
4658
4659       read(inp,*) nset,nfrag,npair,nfrag_back
4660
4661 !el from module energy
4662 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
4663       if(.not.allocated(wfrag_back)) then
4664         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4665         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4666
4667         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4668         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
4669
4670         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4671         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
4672       endif
4673
4674       if(me.eq.king.or..not.out1file) &
4675        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4676         " nfrag_back",nfrag_back
4677       do iset=1,nset
4678          read(inp,*) mset(iset)
4679        do i=1,nfrag
4680          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4681            qinfrag(i,iset)
4682          if(me.eq.king.or..not.out1file) &
4683           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4684            ifrag(2,i,iset), qinfrag(i,iset)
4685        enddo
4686        do i=1,npair
4687         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4688           qinpair(i,iset)
4689         if(me.eq.king.or..not.out1file) &
4690          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4691           ipair(2,i,iset), qinpair(i,iset)
4692        enddo
4693        do i=1,nfrag_back
4694         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4695            wfrag_back(3,i,iset),&
4696            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4697         if(me.eq.king.or..not.out1file) &
4698          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4699          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4700        enddo
4701       enddo
4702       return
4703       end subroutine read_fragments
4704 !-----------------------------------------------------------------------------
4705 ! shift.F       io_csa
4706 !-----------------------------------------------------------------------------
4707       subroutine csa_read
4708
4709       use csa_data
4710 !      implicit real*8 (a-h,o-z)
4711 !      include 'DIMENSIONS'
4712 !      include 'COMMON.CSA'
4713 !      include 'COMMON.BANK'
4714 !      include 'COMMON.IOUNITS'
4715 !el local variables
4716 !     integer :: ntf,ik,iw_pdb
4717 !     real(kind=8) :: var,ene
4718
4719       open(icsa_in,file=csa_in,status="old",err=100)
4720        read(icsa_in,*) nconf
4721        read(icsa_in,*) jstart,jend
4722        read(icsa_in,*) nstmax
4723        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4724        read(icsa_in,*) nran0,nran1,irr
4725        read(icsa_in,*) nseed
4726        read(icsa_in,*) ntotal,cut1,cut2
4727        read(icsa_in,*) estop
4728        read(icsa_in,*) icmax,irestart
4729        read(icsa_in,*) ntbankm,dele,difcut
4730        read(icsa_in,*) iref,rmscut,pnccut
4731        read(icsa_in,*) ndiff
4732       close(icsa_in)
4733
4734       return
4735
4736  100  continue
4737       return
4738       end subroutine csa_read
4739 !-----------------------------------------------------------------------------
4740       subroutine initial_write
4741
4742       use csa_data
4743 !      implicit real*8 (a-h,o-z)
4744 !      include 'DIMENSIONS'
4745 !      include 'COMMON.CSA'
4746 !      include 'COMMON.BANK'
4747 !      include 'COMMON.IOUNITS'
4748 !el local variables
4749 !     integer :: ntf,ik,iw_pdb
4750 !     real(kind=8) :: var,ene
4751
4752       open(icsa_seed,file=csa_seed,status="unknown")
4753        write(icsa_seed,*) "seed"
4754       close(31)
4755 #if defined(AIX) || defined(PGI)
4756        open(icsa_history,file=csa_history,status="unknown",&
4757         position="append")
4758 #else
4759        open(icsa_history,file=csa_history,status="unknown",&
4760         access="append")
4761 #endif
4762        write(icsa_history,*) nconf
4763        write(icsa_history,*) jstart,jend
4764        write(icsa_history,*) nstmax
4765        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4766        write(icsa_history,*) nran0,nran1,irr
4767        write(icsa_history,*) nseed
4768        write(icsa_history,*) ntotal,cut1,cut2
4769        write(icsa_history,*) estop
4770        write(icsa_history,*) icmax,irestart
4771        write(icsa_history,*) ntbankm,dele,difcut
4772        write(icsa_history,*) iref,rmscut,pnccut
4773        write(icsa_history,*) ndiff
4774
4775        write(icsa_history,*)
4776       close(icsa_history)
4777
4778       open(icsa_bank1,file=csa_bank1,status="unknown")
4779        write(icsa_bank1,*) 0
4780       close(icsa_bank1)
4781
4782       return
4783       end subroutine initial_write
4784 !-----------------------------------------------------------------------------
4785       subroutine restart_write
4786
4787       use csa_data
4788 !      implicit real*8 (a-h,o-z)
4789 !      include 'DIMENSIONS'
4790 !      include 'COMMON.IOUNITS'
4791 !      include 'COMMON.CSA'
4792 !      include 'COMMON.BANK'
4793 !el local variables
4794 !     integer :: ntf,ik,iw_pdb
4795 !     real(kind=8) :: var,ene
4796
4797 #if defined(AIX) || defined(PGI)
4798        open(icsa_history,file=csa_history,position="append")
4799 #else
4800        open(icsa_history,file=csa_history,access="append")
4801 #endif
4802        write(icsa_history,*)
4803        write(icsa_history,*) "This is restart"
4804        write(icsa_history,*)
4805        write(icsa_history,*) nconf
4806        write(icsa_history,*) jstart,jend
4807        write(icsa_history,*) nstmax
4808        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4809        write(icsa_history,*) nran0,nran1,irr
4810        write(icsa_history,*) nseed
4811        write(icsa_history,*) ntotal,cut1,cut2
4812        write(icsa_history,*) estop
4813        write(icsa_history,*) icmax,irestart
4814        write(icsa_history,*) ntbankm,dele,difcut
4815        write(icsa_history,*) iref,rmscut,pnccut
4816        write(icsa_history,*) ndiff
4817        write(icsa_history,*)
4818        write(icsa_history,*) "irestart is: ", irestart
4819
4820        write(icsa_history,*)
4821       close(icsa_history)
4822
4823       return
4824       end subroutine restart_write
4825 !-----------------------------------------------------------------------------
4826 ! test.F
4827 !-----------------------------------------------------------------------------
4828       subroutine write_pdb(npdb,titelloc,ee)
4829
4830 !      implicit real*8 (a-h,o-z)
4831 !      include 'DIMENSIONS'
4832 !      include 'COMMON.IOUNITS'
4833       character(len=50) :: titelloc1
4834       character*(*) :: titelloc
4835       character(len=3) :: zahl
4836       character(len=5) :: liczba5
4837       real(kind=8) :: ee
4838       integer :: npdb   !,ilen
4839 !el      external ilen
4840 !el local variables
4841       integer :: lenpre
4842 !     real(kind=8) :: var,ene
4843
4844       titelloc1=titelloc
4845       lenpre=ilen(prefix)
4846       if (npdb.lt.1000) then
4847        call numstr(npdb,zahl)
4848        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4849       else
4850         if (npdb.lt.10000) then
4851          write(liczba5,'(i1,i4)') 0,npdb
4852         else
4853          write(liczba5,'(i5)') npdb
4854         endif
4855         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4856       endif
4857       call pdbout(ee,titelloc1,ipdb)
4858       close(ipdb)
4859       return
4860       end subroutine write_pdb
4861 !-----------------------------------------------------------------------------
4862 ! thread.F
4863 !-----------------------------------------------------------------------------
4864       subroutine write_thread_summary
4865 ! Thread the sequence through a database of known structures
4866       use control_data, only: refstr
4867 !      use geometry
4868       use energy_data, only: n_ene_comp
4869       use compare_data
4870 !      implicit real*8 (a-h,o-z)
4871 !      include 'DIMENSIONS'
4872 #ifdef MPI
4873       use MPI_data      !include 'COMMON.INFO'
4874       include 'mpif.h'
4875 #endif
4876 !      include 'COMMON.CONTROL'
4877 !      include 'COMMON.CHAIN'
4878 !      include 'COMMON.DBASE'
4879 !      include 'COMMON.INTERACT'
4880 !      include 'COMMON.VAR'
4881 !      include 'COMMON.THREAD'
4882 !      include 'COMMON.FFIELD'
4883 !      include 'COMMON.SBRIDGE'
4884 !      include 'COMMON.HEADER'
4885 !      include 'COMMON.NAMES'
4886 !      include 'COMMON.IOUNITS'
4887 !      include 'COMMON.TIME1'
4888
4889       integer,dimension(maxthread) :: ip
4890       real(kind=8),dimension(0:n_ene) :: energia
4891 !el local variables
4892       integer :: i,j,ii,jj,ipj,ik,kk,ist
4893       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4894
4895       write (iout,'(30x,a/)') &
4896        '  *********** Summary threading statistics ************'
4897       write (iout,'(a)') 'Initial energies:'
4898       write (iout,'(a4,2x,a12,14a14,3a8)') &
4899        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4900        'RMSnat','NatCONT','NNCONT','RMS'
4901 ! Energy sort patterns
4902       do i=1,nthread
4903         ip(i)=i
4904       enddo
4905       do i=1,nthread-1
4906         enet=ener(n_ene-1,ip(i))
4907         jj=i
4908         do j=i+1,nthread
4909           if (ener(n_ene-1,ip(j)).lt.enet) then
4910             jj=j
4911             enet=ener(n_ene-1,ip(j))
4912           endif
4913         enddo
4914         if (jj.ne.i) then
4915           ipj=ip(jj)
4916           ip(jj)=ip(i)
4917           ip(i)=ipj
4918         endif
4919       enddo
4920       do ik=1,nthread
4921         i=ip(ik)
4922         ii=ipatt(1,i)
4923         ist=nres_base(2,ii)+ipatt(2,i)
4924         do kk=1,n_ene_comp
4925           energia(i)=ener0(kk,i)
4926         enddo
4927         etot=ener0(n_ene_comp+1,i)
4928         rmsnat=ener0(n_ene_comp+2,i)
4929         rms=ener0(n_ene_comp+3,i)
4930         frac=ener0(n_ene_comp+4,i)
4931         frac_nn=ener0(n_ene_comp+5,i)
4932
4933         if (refstr) then
4934         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4935         i,str_nam(ii),ist+1,&
4936         (energia(print_order(kk)),kk=1,nprint_ene),&
4937         etot,rmsnat,frac,frac_nn,rms
4938         else
4939         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4940         i,str_nam(ii),ist+1,&
4941         (energia(print_order(kk)),kk=1,nprint_ene),etot
4942         endif
4943       enddo
4944       write (iout,'(//a)') 'Final energies:'
4945       write (iout,'(a4,2x,a12,17a14,3a8)') &
4946        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4947        'RMSnat','NatCONT','NNCONT','RMS'
4948       do ik=1,nthread
4949         i=ip(ik)
4950         ii=ipatt(1,i)
4951         ist=nres_base(2,ii)+ipatt(2,i)
4952         do kk=1,n_ene_comp
4953           energia(kk)=ener(kk,ik)
4954         enddo
4955         etot=ener(n_ene_comp+1,i)
4956         rmsnat=ener(n_ene_comp+2,i)
4957         rms=ener(n_ene_comp+3,i)
4958         frac=ener(n_ene_comp+4,i)
4959         frac_nn=ener(n_ene_comp+5,i)
4960         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4961         i,str_nam(ii),ist+1,&
4962         (energia(print_order(kk)),kk=1,nprint_ene),&
4963         etot,rmsnat,frac,frac_nn,rms
4964       enddo
4965       write (iout,'(/a/)') 'IEXAM array:'
4966       write (iout,'(i5)') nexcl
4967       do i=1,nexcl
4968         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4969       enddo
4970       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4971        'Max. time for threading step ',max_time_for_thread,&
4972        'Average time for threading step: ',ave_time_for_thread
4973       return
4974       end subroutine write_thread_summary
4975 !-----------------------------------------------------------------------------
4976       subroutine write_stat_thread(ithread,ipattern,ist)
4977
4978       use energy_data, only: n_ene_comp
4979       use compare_data
4980 !      implicit real*8 (a-h,o-z)
4981 !      include "DIMENSIONS"
4982 !      include "COMMON.CONTROL"
4983 !      include "COMMON.IOUNITS"
4984 !      include "COMMON.THREAD"
4985 !      include "COMMON.FFIELD"
4986 !      include "COMMON.DBASE"
4987 !      include "COMMON.NAMES"
4988       real(kind=8),dimension(0:n_ene) :: energia
4989 !el local variables
4990       integer :: ithread,ipattern,ist,i
4991       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4992
4993 #if defined(AIX) || defined(PGI)
4994       open(istat,file=statname,position='append')
4995 #else
4996       open(istat,file=statname,access='append')
4997 #endif
4998       do i=1,n_ene_comp
4999         energia(i)=ener(i,ithread)
5000       enddo
5001       etot=ener(n_ene_comp+1,ithread)
5002       rmsnat=ener(n_ene_comp+2,ithread)
5003       rms=ener(n_ene_comp+3,ithread)
5004       frac=ener(n_ene_comp+4,ithread)
5005       frac_nn=ener(n_ene_comp+5,ithread)
5006       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5007         ithread,str_nam(ipattern),ist+1,&
5008         (energia(print_order(i)),i=1,nprint_ene),&
5009         etot,rmsnat,frac,frac_nn,rms
5010       close (istat)
5011       return
5012       end subroutine write_stat_thread
5013 !-----------------------------------------------------------------------------
5014 #endif
5015 !-----------------------------------------------------------------------------
5016       end module io_config