source/unres/io_config.f90

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