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