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(ntyp)) !(ntyp)
 780       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 781       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 782       allocate(msc(ntyp+1)) !(ntyp+1)
 783       allocate(isc(ntyp+1)) !(ntyp+1)
 784       allocate(restok(ntyp+1)) !(ntyp+1)
 785       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 786       allocate(long_r_sidechain(ntyp))
 787       allocate(short_r_sidechain(ntyp))
 788       dsc(:)=0.0d0
 789       dsc_inv(:)=0.0d0
 790
 791 #ifdef CRYST_BOND
 792       read (ibond,*) vbldp0,akp,mp,ip,pstok
 793       do i=1,ntyp
 794         nbondterm(i)=1
 795         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
 796         dsc(i) = vbldsc0(1,i)
 797         if (i.eq.10) then
 798           dsc_inv(i)=0.0D0
 799         else
 800           dsc_inv(i)=1.0D0/dsc(i)
 801         endif
 802       enddo
 803 #else
 804       read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
 805       do i=1,ntyp
 806         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
 807          j=1,nbondterm(i)),msc(i),isc(i),restok(i)
 808         dsc(i) = vbldsc0(1,i)
 809         if (i.eq.10) then
 810           dsc_inv(i)=0.0D0
 811         else
 812           dsc_inv(i)=1.0D0/dsc(i)
 813         endif
 814       enddo
 815 #endif
 816       if (lprint) then
 817         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
 818         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
 819          'inertia','Pstok'
 820         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
 821         do i=1,ntyp
 822           write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
 823             vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
 824           do j=2,nbondterm(i)
 825             write (iout,'(13x,3f10.5)') &
 826               vbldsc0(j,i),aksc(j,i),abond0(j,i)
 827           enddo
 828         enddo
 829       endif
 830 !----------------------------------------------------
 831       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
 832       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))      !(-ntyp:ntyp)
 833       allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
 834       allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
 835       allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
 836       allocate(gthet(3,-ntyp:ntyp))     !(3,-ntyp:ntyp)
 837
 838       a0thet(:)=0.0D0
 839       athet(:,:,:,:)=0.0D0
 840       bthet(:,:,:,:)=0.0D0
 841       polthet(:,:)=0.0D0
 842       gthet(:,:)=0.0D0
 843       theta0(:)=0.0D0
 844       sig0(:)=0.0D0
 845       sigc0(:)=0.0D0
 846       allocate(liptranene(ntyp))
 847 !C reading lipid parameters
 848       write (iout,*) "iliptranpar",iliptranpar
 849       call flush(iout)
 850        read(iliptranpar,*) pepliptran
 851        print *,pepliptran
 852        do i=1,ntyp
 853        read(iliptranpar,*) liptranene(i)
 854         print *,liptranene(i)
 855        enddo
 856        close(iliptranpar)
 857
 858 #ifdef CRYST_THETA
 859 !
 860 ! Read the parameters of the probability distribution/energy expression
 861 ! of the virtual-bond valence angles theta
 862 !
 863       do i=1,ntyp
 864         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
 865           (bthet(j,i,1,1),j=1,2)
 866         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 867         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 868         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 869         sigc0(i)=sigc0(i)**2
 870       enddo
 871       do i=1,ntyp
 872       athet(1,i,1,-1)=athet(1,i,1,1)
 873       athet(2,i,1,-1)=athet(2,i,1,1)
 874       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 875       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 876       athet(1,i,-1,1)=-athet(1,i,1,1)
 877       athet(2,i,-1,1)=-athet(2,i,1,1)
 878       bthet(1,i,-1,1)=bthet(1,i,1,1)
 879       bthet(2,i,-1,1)=bthet(2,i,1,1)
 880       enddo
 881       do i=-ntyp,-1
 882       a0thet(i)=a0thet(-i)
 883       athet(1,i,-1,-1)=athet(1,-i,1,1)
 884       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 885       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 886       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 887       athet(1,i,-1,1)=athet(1,-i,1,1)
 888       athet(2,i,-1,1)=-athet(2,-i,1,1)
 889       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 890       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 891       athet(1,i,1,-1)=-athet(1,-i,1,1)
 892       athet(2,i,1,-1)=athet(2,-i,1,1)
 893       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 894       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 895       theta0(i)=theta0(-i)
 896       sig0(i)=sig0(-i)
 897       sigc0(i)=sigc0(-i)
 898        do j=0,3
 899         polthet(j,i)=polthet(j,-i)
 900        enddo
 901        do j=1,3
 902          gthet(j,i)=gthet(j,-i)
 903        enddo
 904       enddo
 905
 906       close (ithep)
 907       if (lprint) then
 908       if (.not.LaTeX) then
 909         write (iout,'(a)') &
 910           'Parameters of the virtual-bond valence angles:'
 911         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
 912        '    ATHETA0   ','         A1   ','        A2    ',&
 913        '        B1    ','         B2   '
 914         do i=1,ntyp
 915           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
 916               a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
 917         enddo
 918         write (iout,'(/a/9x,5a/79(1h-))') &
 919        'Parameters of the expression for sigma(theta_c):',&
 920        '     ALPH0    ','      ALPH1   ','     ALPH2    ',&
 921        '     ALPH3    ','    SIGMA0C   '
 922         do i=1,ntyp
 923           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
 924             (polthet(j,i),j=0,3),sigc0(i)
 925         enddo
 926         write (iout,'(/a/9x,5a/79(1h-))') &
 927        'Parameters of the second gaussian:',&
 928        '    THETA0    ','     SIGMA0   ','        G1    ',&
 929        '        G2    ','         G3   '
 930         do i=1,ntyp
 931           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
 932              sig0(i),(gthet(j,i),j=1,3)
 933         enddo
 934        else
 935         write (iout,'(a)') &
 936           'Parameters of the virtual-bond valence angles:'
 937         write (iout,'(/a/9x,5a/79(1h-))') &
 938        'Coefficients of expansion',&
 939        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
 940        '   b1*10^1    ','    b2*10^1   '
 941         do i=1,ntyp
 942           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
 943               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
 944               (10*bthet(j,i,1,1),j=1,2)
 945         enddo
 946         write (iout,'(/a/9x,5a/79(1h-))') &
 947        'Parameters of the expression for sigma(theta_c):',&
 948        ' alpha0       ','  alph1       ',' alph2        ',&
 949        ' alhp3        ','   sigma0c    '
 950         do i=1,ntyp
 951           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
 952             (polthet(j,i),j=0,3),sigc0(i)
 953         enddo
 954         write (iout,'(/a/9x,5a/79(1h-))') &
 955        'Parameters of the second gaussian:',&
 956        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
 957        '        G2    ','   G3*10^1    '
 958         do i=1,ntyp
 959           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
 960              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 961         enddo
 962       endif
 963       endif
 964 #else
 965 !
 966 ! Read the parameters of Utheta determined from ab initio surfaces
 967 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 968 !
 969       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
 970         ntheterm3,nsingle,ndouble
 971       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 972
 973 !----------------------------------------------------
 974       allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
 975       allocate(aa0thet(-nthetyp-1:nthetyp+1,&
 976         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 977 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 978       allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
 979         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 980 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 981 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 982       allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 983         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 984       allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 985         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 986       allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 987         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 988       allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
 989         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 990 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
 991 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 992       allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 993         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 994       allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
 995         -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
 996 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
 997 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 998
 999       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1000       do i=-ntyp1,-1
1001         ithetyp(i)=-ithetyp(-i)
1002       enddo
1003
1004       aa0thet(:,:,:,:)=0.0d0
1005       aathet(:,:,:,:,:)=0.0d0
1006       bbthet(:,:,:,:,:,:)=0.0d0
1007       ccthet(:,:,:,:,:,:)=0.0d0
1008       ddthet(:,:,:,:,:,:)=0.0d0
1009       eethet(:,:,:,:,:,:)=0.0d0
1010       ffthet(:,:,:,:,:,:,:)=0.0d0
1011       ggthet(:,:,:,:,:,:,:)=0.0d0
1012
1013 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1014       do iblock=1,2
1015 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1016 ! VAR:1=non-glicyne non-proline 2=proline
1017 ! VAR:negative values for D-aminoacid
1018       do i=0,nthetyp
1019         do j=-nthetyp,nthetyp
1020           do k=-nthetyp,nthetyp
1021             read (ithep,'(6a)',end=111,err=111) res1
1022             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1023 ! VAR: aa0thet is variable describing the average value of Foureir
1024 ! VAR: expansion series
1025 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1026 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1027 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1028             read (ithep,*,end=111,err=111) &
1029               (aathet(l,i,j,k,iblock),l=1,ntheterm)
1030             read (ithep,*,end=111,err=111) &
1031              ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1032               (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1033               (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1034               (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1035               ll=1,ntheterm2)
1036             read (ithep,*,end=111,err=111) &
1037             (((ffthet(llll,lll,ll,i,j,k,iblock),&
1038                ffthet(lll,llll,ll,i,j,k,iblock),&
1039                ggthet(llll,lll,ll,i,j,k,iblock),&
1040                ggthet(lll,llll,ll,i,j,k,iblock),&
1041                llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1042           enddo
1043         enddo
1044       enddo
1045 !
1046 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1047 ! coefficients of theta-and-gamma-dependent terms are zero.
1048 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1049 ! RECOMENTDED AFTER VERSION 3.3)
1050 !      do i=1,nthetyp
1051 !        do j=1,nthetyp
1052 !          do l=1,ntheterm
1053 !            aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1054 !            aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1055 !          enddo
1056 !          aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1057 !          aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1058 !        enddo
1059 !        do l=1,ntheterm
1060 !          aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1061 !        enddo
1062 !        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1063 !      enddo
1064 !      enddo
1065 ! AND COMMENT THE LOOPS BELOW
1066       do i=1,nthetyp
1067         do j=1,nthetyp
1068           do l=1,ntheterm
1069             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1070             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1071           enddo
1072           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1073           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1074         enddo
1075         do l=1,ntheterm
1076           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1077         enddo
1078         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1079       enddo
1080       enddo       !iblock
1081
1082 ! TILL HERE
1083 ! Substitution for D aminoacids from symmetry.
1084       do iblock=1,2
1085       do i=-nthetyp,0
1086         do j=-nthetyp,nthetyp
1087           do k=-nthetyp,nthetyp
1088            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1089            do l=1,ntheterm
1090            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1091            enddo
1092            do ll=1,ntheterm2
1093             do lll=1,nsingle
1094             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1095             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1096             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1097             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1098             enddo
1099           enddo
1100           do ll=1,ntheterm3
1101            do lll=2,ndouble
1102             do llll=1,lll-1
1103             ffthet(llll,lll,ll,i,j,k,iblock)= &
1104             ffthet(llll,lll,ll,-i,-j,-k,iblock)
1105             ffthet(lll,llll,ll,i,j,k,iblock)= &
1106             ffthet(lll,llll,ll,-i,-j,-k,iblock)
1107             ggthet(llll,lll,ll,i,j,k,iblock)= &
1108             -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1109             ggthet(lll,llll,ll,i,j,k,iblock)= &
1110             -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1111             enddo !ll
1112            enddo  !lll
1113           enddo   !llll
1114          enddo    !k
1115         enddo     !j
1116        enddo      !i
1117       enddo       !iblock
1118 !
1119 ! Control printout of the coefficients of virtual-bond-angle potentials
1120 !
1121       if (lprint) then
1122         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1123         do iblock=1,2
1124         do i=1,nthetyp+1
1125           do j=1,nthetyp+1
1126             do k=1,nthetyp+1
1127               write (iout,'(//4a)') &
1128                'Type ',onelett(i),onelett(j),onelett(k)
1129               write (iout,'(//a,10x,a)') " l","a[l]"
1130               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1131               write (iout,'(i2,1pe15.5)') &
1132                  (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1133             do l=1,ntheterm2
1134               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1135                 "b",l,"c",l,"d",l,"e",l
1136               do m=1,nsingle
1137                 write (iout,'(i2,4(1pe15.5))') m,&
1138                 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1139                 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1140               enddo
1141             enddo
1142             do l=1,ntheterm3
1143               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1144                 "f+",l,"f-",l,"g+",l,"g-",l
1145               do m=2,ndouble
1146                 do n=1,m-1
1147                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1148                     ffthet(n,m,l,i,j,k,iblock),&
1149                     ffthet(m,n,l,i,j,k,iblock),&
1150                     ggthet(n,m,l,i,j,k,iblock),&
1151                     ggthet(m,n,l,i,j,k,iblock)
1152                 enddo   !n
1153               enddo     !m
1154             enddo       !l
1155           enddo         !k
1156         enddo           !j
1157       enddo             !i
1158       enddo
1159       call flush(iout)
1160       endif
1161       write (2,*) "Start reading THETA_PDB",ithep_pdb
1162       do i=1,ntyp
1163 !      write (2,*) 'i=',i
1164         read (ithep_pdb,*,err=111,end=111) &
1165            a0thet(i),(athet(j,i,1,1),j=1,2),&
1166           (bthet(j,i,1,1),j=1,2)
1167         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1168         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1169         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1170         sigc0(i)=sigc0(i)**2
1171       enddo
1172       do i=1,ntyp
1173       athet(1,i,1,-1)=athet(1,i,1,1)
1174       athet(2,i,1,-1)=athet(2,i,1,1)
1175       bthet(1,i,1,-1)=-bthet(1,i,1,1)
1176       bthet(2,i,1,-1)=-bthet(2,i,1,1)
1177       athet(1,i,-1,1)=-athet(1,i,1,1)
1178       athet(2,i,-1,1)=-athet(2,i,1,1)
1179       bthet(1,i,-1,1)=bthet(1,i,1,1)
1180       bthet(2,i,-1,1)=bthet(2,i,1,1)
1181       enddo
1182       do i=-ntyp,-1
1183       a0thet(i)=a0thet(-i)
1184       athet(1,i,-1,-1)=athet(1,-i,1,1)
1185       athet(2,i,-1,-1)=-athet(2,-i,1,1)
1186       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1187       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1188       athet(1,i,-1,1)=athet(1,-i,1,1)
1189       athet(2,i,-1,1)=-athet(2,-i,1,1)
1190       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1191       bthet(2,i,-1,1)=bthet(2,-i,1,1)
1192       athet(1,i,1,-1)=-athet(1,-i,1,1)
1193       athet(2,i,1,-1)=athet(2,-i,1,1)
1194       bthet(1,i,1,-1)=bthet(1,-i,1,1)
1195       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1196       theta0(i)=theta0(-i)
1197       sig0(i)=sig0(-i)
1198       sigc0(i)=sigc0(-i)
1199        do j=0,3
1200         polthet(j,i)=polthet(j,-i)
1201        enddo
1202        do j=1,3
1203          gthet(j,i)=gthet(j,-i)
1204        enddo
1205       enddo
1206       write (2,*) "End reading THETA_PDB"
1207       close (ithep_pdb)
1208 #endif
1209       close(ithep)
1210
1211 !-------------------------------------------
1212       allocate(nlob(ntyp1)) !(ntyp1)
1213       allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1214       allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1215       allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1216
1217       bsc(:,:)=0.0D0
1218       nlob(:)=0
1219       nlob(:)=0
1220       dsc(:)=0.0D0
1221       censc(:,:,:)=0.0D0
1222       gaussc(:,:,:,:)=0.0D0
1223
1224 #ifdef CRYST_SC
1225 !
1226 ! Read the parameters of the probability distribution/energy expression
1227 ! of the side chains.
1228 !
1229       do i=1,ntyp
1230         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1231         if (i.eq.10) then
1232           dsc_inv(i)=0.0D0
1233         else
1234           dsc_inv(i)=1.0D0/dsc(i)
1235         endif
1236         if (i.ne.10) then
1237         do j=1,nlob(i)
1238           do k=1,3
1239             do l=1,3
1240               blower(l,k,j)=0.0D0
1241             enddo
1242           enddo
1243         enddo
1244         bsc(1,i)=0.0D0
1245         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1246           ((blower(k,l,1),l=1,k),k=1,3)
1247         censc(1,1,-i)=censc(1,1,i)
1248         censc(2,1,-i)=censc(2,1,i)
1249         censc(3,1,-i)=-censc(3,1,i)
1250         do j=2,nlob(i)
1251           read (irotam,*,end=112,err=112) bsc(j,i)
1252           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1253                                        ((blower(k,l,j),l=1,k),k=1,3)
1254         censc(1,j,-i)=censc(1,j,i)
1255         censc(2,j,-i)=censc(2,j,i)
1256         censc(3,j,-i)=-censc(3,j,i)
1257 ! BSC is amplitude of Gaussian
1258         enddo
1259         do j=1,nlob(i)
1260           do k=1,3
1261             do l=1,k
1262               akl=0.0D0
1263               do m=1,3
1264                 akl=akl+blower(k,m,j)*blower(l,m,j)
1265               enddo
1266               gaussc(k,l,j,i)=akl
1267               gaussc(l,k,j,i)=akl
1268              if (((k.eq.3).and.(l.ne.3)) &
1269               .or.((l.eq.3).and.(k.ne.3))) then
1270                 gaussc(k,l,j,-i)=-akl
1271                 gaussc(l,k,j,-i)=-akl
1272               else
1273                 gaussc(k,l,j,-i)=akl
1274                 gaussc(l,k,j,-i)=akl
1275               endif
1276             enddo
1277           enddo
1278         enddo
1279         endif
1280       enddo
1281       close (irotam)
1282       if (lprint) then
1283         write (iout,'(/a)') 'Parameters of side-chain local geometry'
1284         do i=1,ntyp
1285           nlobi=nlob(i)
1286           if (nlobi.gt.0) then
1287             if (LaTeX) then
1288               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1289                ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1290                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1291                                    'log h',(bsc(j,i),j=1,nlobi)
1292                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1293               'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1294               do k=1,3
1295                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1296                        ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1297               enddo
1298             else
1299               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1300               write (iout,'(a,f10.4,4(16x,f10.4))') &
1301                                    'Center  ',(bsc(j,i),j=1,nlobi)
1302               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1303                  j=1,nlobi)
1304               write (iout,'(a)')
1305             endif
1306           endif
1307         enddo
1308       endif
1309 #else
1310 !
1311 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1312 ! added by Urszula Kozlowska 07/11/2007
1313 !
1314 !el Maximum number of SC local term fitting function coefficiants
1315 !el      integer,parameter :: maxsccoef=65
1316
1317       allocate(sc_parmin(65,ntyp))      !(maxsccoef,ntyp)
1318
1319       do i=1,ntyp
1320         read (irotam,*,end=112,err=112)
1321        if (i.eq.10) then
1322          read (irotam,*,end=112,err=112)
1323        else
1324          do j=1,65
1325            read(irotam,*,end=112,err=112) sc_parmin(j,i)
1326          enddo
1327        endif
1328       enddo
1329 !
1330 ! Read the parameters of the probability distribution/energy expression
1331 ! of the side chains.
1332 !
1333       write (2,*) "Start reading ROTAM_PDB"
1334       do i=1,ntyp
1335         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1336         if (i.eq.10) then
1337           dsc_inv(i)=0.0D0
1338         else
1339           dsc_inv(i)=1.0D0/dsc(i)
1340         endif
1341         if (i.ne.10) then
1342         do j=1,nlob(i)
1343           do k=1,3
1344             do l=1,3
1345               blower(l,k,j)=0.0D0
1346             enddo
1347           enddo
1348         enddo
1349         bsc(1,i)=0.0D0
1350         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1351           ((blower(k,l,1),l=1,k),k=1,3)
1352         do j=2,nlob(i)
1353           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1354           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1355                                        ((blower(k,l,j),l=1,k),k=1,3)
1356         enddo
1357         do j=1,nlob(i)
1358           do k=1,3
1359             do l=1,k
1360               akl=0.0D0
1361               do m=1,3
1362                 akl=akl+blower(k,m,j)*blower(l,m,j)
1363               enddo
1364               gaussc(k,l,j,i)=akl
1365               gaussc(l,k,j,i)=akl
1366             enddo
1367           enddo
1368         enddo
1369         endif
1370       enddo
1371       close (irotam_pdb)
1372       write (2,*) "End reading ROTAM_PDB"
1373 #endif
1374       close(irotam)
1375
1376 #ifdef CRYST_TOR
1377 !
1378 ! Read torsional parameters in old format
1379 !
1380       allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1381
1382       read (itorp,*,end=113,err=113) ntortyp,nterm_old
1383       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1384       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1385
1386 !el from energy module--------
1387       allocate(v1(nterm_old,ntortyp,ntortyp))
1388       allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1389 !el---------------------------
1390       do i=1,ntortyp
1391         do j=1,ntortyp
1392           read (itorp,'(a)')
1393           do k=1,nterm_old
1394             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1395           enddo
1396         enddo
1397       enddo
1398       close (itorp)
1399       if (lprint) then
1400         write (iout,'(/a/)') 'Torsional constants:'
1401         do i=1,ntortyp
1402           do j=1,ntortyp
1403             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1404             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1405           enddo
1406         enddo
1407       endif
1408 #else
1409 !
1410 ! Read torsional parameters
1411 !
1412       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1413       read (itorp,*,end=113,err=113) ntortyp
1414 !el from energy module---------
1415       allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1416       allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1417
1418       allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1419       allocate(vlor2(maxlor,ntortyp,ntortyp))
1420       allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1421       allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1422
1423       allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1424       allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1425 !el---------------------------
1426       nterm(:,:,:)=0
1427       nlor(:,:,:)=0
1428 !el---------------------------
1429
1430       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1431       do i=-ntyp,-1
1432        itortyp(i)=-itortyp(-i)
1433       enddo
1434       itortyp(ntyp1)=ntortyp
1435       itortyp(-ntyp1)=-ntortyp
1436       do iblock=1,2
1437       write (iout,*) 'ntortyp',ntortyp
1438       do i=0,ntortyp-1
1439         do j=-ntortyp+1,ntortyp-1
1440           read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1441                 nlor(i,j,iblock)
1442           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1443           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1444           v0ij=0.0d0
1445           si=-1.0d0
1446           do k=1,nterm(i,j,iblock)
1447             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1448             v2(k,i,j,iblock)
1449             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1450             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1451             v0ij=v0ij+si*v1(k,i,j,iblock)
1452             si=-si
1453 !         write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1454 !         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1455 !      v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1456           enddo
1457           do k=1,nlor(i,j,iblock)
1458             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1459               vlor2(k,i,j),vlor3(k,i,j)
1460             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1461           enddo
1462           v0(i,j,iblock)=v0ij
1463           v0(-i,-j,iblock)=v0ij
1464         enddo
1465       enddo
1466       enddo
1467       close (itorp)
1468       if (lprint) then
1469         write (iout,'(/a/)') 'Torsional constants:'
1470         do iblock=1,2
1471         do i=-ntortyp,ntortyp
1472           do j=-ntortyp,ntortyp
1473             write (iout,*) 'ityp',i,' jtyp',j
1474             write (iout,*) 'Fourier constants'
1475             do k=1,nterm(i,j,iblock)
1476               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1477               v2(k,i,j,iblock)
1478             enddo
1479             write (iout,*) 'Lorenz constants'
1480             do k=1,nlor(i,j,iblock)
1481               write (iout,'(3(1pe15.5))') &
1482                vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1483             enddo
1484           enddo
1485         enddo
1486         enddo
1487       endif
1488 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1489 !
1490 ! 6/23/01 Read parameters for double torsionals
1491 !
1492 !el from energy module------------
1493       allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1494       allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1495 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1496       allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1497       allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1498         !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1499       allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1500       allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1501         !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1502 !---------------------------------
1503
1504       do iblock=1,2
1505       do i=0,ntortyp-1
1506         do j=-ntortyp+1,ntortyp-1
1507           do k=-ntortyp+1,ntortyp-1
1508             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1509 !              write (iout,*) "OK onelett",
1510 !     &         i,j,k,t1,t2,t3
1511
1512             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1513               .or. t3.ne.toronelet(k)) then
1514              write (iout,*) "Error in double torsional parameter file",&
1515                i,j,k,t1,t2,t3
1516 #ifdef MPI
1517               call MPI_Finalize(Ierror)
1518 #endif
1519                stop "Error in double torsional parameter file"
1520             endif
1521            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1522                ntermd_2(i,j,k,iblock)
1523             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1524             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1525             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1526                ntermd_1(i,j,k,iblock))
1527             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1528                ntermd_1(i,j,k,iblock))
1529             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1530                ntermd_1(i,j,k,iblock))
1531             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1532                ntermd_1(i,j,k,iblock))
1533 ! Martix of D parameters for one dimesional foureir series
1534             do l=1,ntermd_1(i,j,k,iblock)
1535              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1536              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1537              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1538              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1539 !            write(iout,*) "whcodze" ,
1540 !     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1541             enddo
1542             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1543                v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1544                v2s(m,l,i,j,k,iblock),&
1545                m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1546 ! Martix of D parameters for two dimesional fourier series
1547             do l=1,ntermd_2(i,j,k,iblock)
1548              do m=1,l-1
1549              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1550              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1551              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1552              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1553              enddo!m
1554             enddo!l
1555           enddo!k
1556         enddo!j
1557       enddo!i
1558       enddo!iblock
1559       if (lprint) then
1560       write (iout,*)
1561       write (iout,*) 'Constants for double torsionals'
1562       do iblock=1,2
1563       do i=0,ntortyp-1
1564         do j=-ntortyp+1,ntortyp-1
1565           do k=-ntortyp+1,ntortyp-1
1566             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1567               ' nsingle',ntermd_1(i,j,k,iblock),&
1568               ' ndouble',ntermd_2(i,j,k,iblock)
1569             write (iout,*)
1570             write (iout,*) 'Single angles:'
1571             do l=1,ntermd_1(i,j,k,iblock)
1572               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1573                  v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1574                  v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1575                  v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1576             enddo
1577             write (iout,*)
1578             write (iout,*) 'Pairs of angles:'
1579             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1580             do l=1,ntermd_2(i,j,k,iblock)
1581               write (iout,'(i5,20f10.5)') &
1582                l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1583             enddo
1584             write (iout,*)
1585             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1586             do l=1,ntermd_2(i,j,k,iblock)
1587               write (iout,'(i5,20f10.5)') &
1588                l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1589                (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1590             enddo
1591             write (iout,*)
1592           enddo
1593         enddo
1594       enddo
1595       enddo
1596       endif
1597 #endif
1598 ! Read of Side-chain backbone correlation parameters
1599 ! Modified 11 May 2012 by Adasko
1600 !CC
1601 !
1602       read (isccor,*,end=119,err=119) nsccortyp
1603
1604 !el from module energy-------------
1605       allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1606       allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1607       allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1608       allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp))   !(maxterm_sccor,20,20)
1609 !-----------------------------------
1610 #ifdef SCCORPDB
1611 !el from module energy-------------
1612       allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1613
1614       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1615       do i=-ntyp,-1
1616         isccortyp(i)=-isccortyp(-i)
1617       enddo
1618       iscprol=isccortyp(20)
1619 !      write (iout,*) 'ntortyp',ntortyp
1620       maxinter=3
1621 !c maxinter is maximum interaction sites
1622 !el from module energy---------
1623       allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1624       allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1625                -nsccortyp:nsccortyp))
1626       allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1627                -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1628       allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1629                -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1630 !-----------------------------------
1631       do l=1,maxinter
1632       do i=1,nsccortyp
1633         do j=1,nsccortyp
1634           read (isccor,*,end=119,err=119) &
1635       nterm_sccor(i,j),nlor_sccor(i,j)
1636           v0ijsccor=0.0d0
1637           v0ijsccor1=0.0d0
1638           v0ijsccor2=0.0d0
1639           v0ijsccor3=0.0d0
1640           si=-1.0d0
1641           nterm_sccor(-i,j)=nterm_sccor(i,j)
1642           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1643           nterm_sccor(i,-j)=nterm_sccor(i,j)
1644           do k=1,nterm_sccor(i,j)
1645             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1646            v2sccor(k,l,i,j)
1647             if (j.eq.iscprol) then
1648              if (i.eq.isccortyp(10)) then
1649              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1650              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1651              else
1652              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1653                               +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1654              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1655                               +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1656              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1657              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1658              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1659              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1660              endif
1661             else
1662              if (i.eq.isccortyp(10)) then
1663              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1664              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1665              else
1666                if (j.eq.isccortyp(10)) then
1667              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1668              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1669                else
1670              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1671              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1672              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1673              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1674              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1675              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1676                 endif
1677                endif
1678             endif
1679             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1680             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1681             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1682             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1683             si=-si
1684           enddo
1685           do k=1,nlor_sccor(i,j)
1686             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1687               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1688             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1689       (1+vlor3sccor(k,i,j)**2)
1690           enddo
1691           v0sccor(l,i,j)=v0ijsccor
1692           v0sccor(l,-i,j)=v0ijsccor1
1693           v0sccor(l,i,-j)=v0ijsccor2
1694           v0sccor(l,-i,-j)=v0ijsccor3
1695         enddo
1696       enddo
1697       enddo
1698       close (isccor)
1699 #else
1700 !el from module energy-------------
1701       allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1702
1703       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1704 !      write (iout,*) 'ntortyp',ntortyp
1705       maxinter=3
1706 !c maxinter is maximum interaction sites
1707 !el from module energy---------
1708       allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1709       allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1710       allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1711       allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1712 !-----------------------------------
1713       do l=1,maxinter
1714       do i=1,nsccortyp
1715         do j=1,nsccortyp
1716           read (isccor,*,end=119,err=119) &
1717        nterm_sccor(i,j),nlor_sccor(i,j)
1718           v0ijsccor=0.0d0
1719           si=-1.0d0
1720
1721           do k=1,nterm_sccor(i,j)
1722             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1723            v2sccor(k,l,i,j)
1724             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1725             si=-si
1726           enddo
1727           do k=1,nlor_sccor(i,j)
1728             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1729               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1730             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1731       (1+vlor3sccor(k,i,j)**2)
1732           enddo
1733           v0sccor(l,i,j)=v0ijsccor !el ,iblock
1734         enddo
1735       enddo
1736       enddo
1737       close (isccor)
1738
1739 #endif
1740       if (lprint) then
1741         write (iout,'(/a/)') 'Torsional constants:'
1742         do i=1,nsccortyp
1743           do j=1,nsccortyp
1744             write (iout,*) 'ityp',i,' jtyp',j
1745             write (iout,*) 'Fourier constants'
1746             do k=1,nterm_sccor(i,j)
1747       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1748             enddo
1749             write (iout,*) 'Lorenz constants'
1750             do k=1,nlor_sccor(i,j)
1751               write (iout,'(3(1pe15.5))') &
1752                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1753             enddo
1754           enddo
1755         enddo
1756       endif
1757
1758 !
1759 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1760 !         interaction energy of the Gly, Ala, and Pro prototypes.
1761 !
1762
1763       if (lprint) then
1764         write (iout,*)
1765         write (iout,*) "Coefficients of the cumulants"
1766       endif
1767       read (ifourier,*) nloctyp
1768 !write(iout,*) "nloctyp",nloctyp
1769 !el from module energy-------
1770       allocate(b1(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1771       allocate(b2(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1772       allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1773       allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1774       allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1775       allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1776       allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1777       allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1778 ! el
1779         b1(1,:)=0.0d0
1780         b1(2,:)=0.0d0
1781 !--------------------------------
1782
1783       do i=0,nloctyp-1
1784         read (ifourier,*,end=115,err=115)
1785         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1786         if (lprint) then
1787           write (iout,*) 'Type',i
1788           write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1789         endif
1790         B1(1,i)  = b(3)
1791         B1(2,i)  = b(5)
1792         B1(1,-i) = b(3)
1793         B1(2,-i) = -b(5)
1794 !        b1(1,i)=0.0d0
1795 !        b1(2,i)=0.0d0
1796         B1tilde(1,i) = b(3)
1797         B1tilde(2,i) =-b(5)
1798         B1tilde(1,-i) =-b(3)
1799         B1tilde(2,-i) =b(5)
1800 !        b1tilde(1,i)=0.0d0
1801 !        b1tilde(2,i)=0.0d0
1802         B2(1,i)  = b(2)
1803         B2(2,i)  = b(4)
1804         B2(1,-i)  =b(2)
1805         B2(2,-i)  =-b(4)
1806
1807 !        b2(1,i)=0.0d0
1808 !        b2(2,i)=0.0d0
1809         CC(1,1,i)= b(7)
1810         CC(2,2,i)=-b(7)
1811         CC(2,1,i)= b(9)
1812         CC(1,2,i)= b(9)
1813         CC(1,1,-i)= b(7)
1814         CC(2,2,-i)=-b(7)
1815         CC(2,1,-i)=-b(9)
1816         CC(1,2,-i)=-b(9)
1817 !        CC(1,1,i)=0.0d0
1818 !        CC(2,2,i)=0.0d0
1819 !        CC(2,1,i)=0.0d0
1820 !        CC(1,2,i)=0.0d0
1821         Ctilde(1,1,i)=b(7)
1822         Ctilde(1,2,i)=b(9)
1823         Ctilde(2,1,i)=-b(9)
1824         Ctilde(2,2,i)=b(7)
1825         Ctilde(1,1,-i)=b(7)
1826         Ctilde(1,2,-i)=-b(9)
1827         Ctilde(2,1,-i)=b(9)
1828         Ctilde(2,2,-i)=b(7)
1829
1830 !        Ctilde(1,1,i)=0.0d0
1831 !        Ctilde(1,2,i)=0.0d0
1832 !        Ctilde(2,1,i)=0.0d0
1833 !        Ctilde(2,2,i)=0.0d0
1834         DD(1,1,i)= b(6)
1835         DD(2,2,i)=-b(6)
1836         DD(2,1,i)= b(8)
1837         DD(1,2,i)= b(8)
1838         DD(1,1,-i)= b(6)
1839         DD(2,2,-i)=-b(6)
1840         DD(2,1,-i)=-b(8)
1841         DD(1,2,-i)=-b(8)
1842 !        DD(1,1,i)=0.0d0
1843 !        DD(2,2,i)=0.0d0
1844 !        DD(2,1,i)=0.0d0
1845 !        DD(1,2,i)=0.0d0
1846         Dtilde(1,1,i)=b(6)
1847         Dtilde(1,2,i)=b(8)
1848         Dtilde(2,1,i)=-b(8)
1849         Dtilde(2,2,i)=b(6)
1850         Dtilde(1,1,-i)=b(6)
1851         Dtilde(1,2,-i)=-b(8)
1852         Dtilde(2,1,-i)=b(8)
1853         Dtilde(2,2,-i)=b(6)
1854
1855 !        Dtilde(1,1,i)=0.0d0
1856 !        Dtilde(1,2,i)=0.0d0
1857 !        Dtilde(2,1,i)=0.0d0
1858 !        Dtilde(2,2,i)=0.0d0
1859         EE(1,1,i)= b(10)+b(11)
1860         EE(2,2,i)=-b(10)+b(11)
1861         EE(2,1,i)= b(12)-b(13)
1862         EE(1,2,i)= b(12)+b(13)
1863         EE(1,1,-i)= b(10)+b(11)
1864         EE(2,2,-i)=-b(10)+b(11)
1865         EE(2,1,-i)=-b(12)+b(13)
1866         EE(1,2,-i)=-b(12)-b(13)
1867
1868 !        ee(1,1,i)=1.0d0
1869 !        ee(2,2,i)=1.0d0
1870 !        ee(2,1,i)=0.0d0
1871 !        ee(1,2,i)=0.0d0
1872 !        ee(2,1,i)=ee(1,2,i)
1873       enddo
1874       if (lprint) then
1875       do i=1,nloctyp
1876         write (iout,*) 'Type',i
1877         write (iout,*) 'B1'
1878         write(iout,*) B1(1,i),B1(2,i)
1879         write (iout,*) 'B2'
1880         write(iout,*) B2(1,i),B2(2,i)
1881         write (iout,*) 'CC'
1882         do j=1,2
1883           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1884         enddo
1885         write(iout,*) 'DD'
1886         do j=1,2
1887           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1888         enddo
1889         write(iout,*) 'EE'
1890         do j=1,2
1891           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1892         enddo
1893       enddo
1894       endif
1895 !
1896 ! Read electrostatic-interaction parameters
1897 !
1898
1899       if (lprint) then
1900         write (iout,*)
1901         write (iout,'(/a)') 'Electrostatic interaction constants:'
1902         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1903                   'IT','JT','APP','BPP','AEL6','AEL3'
1904       endif
1905       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1906       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1907       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1908       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1909       close (ielep)
1910       do i=1,2
1911         do j=1,2
1912         rri=rpp(i,j)**6
1913         app (i,j)=epp(i,j)*rri*rri
1914         bpp (i,j)=-2.0D0*epp(i,j)*rri
1915         ael6(i,j)=elpp6(i,j)*4.2D0**6
1916         ael3(i,j)=elpp3(i,j)*4.2D0**3
1917 !        lprint=.true.
1918         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1919                           ael6(i,j),ael3(i,j)
1920 !        lprint=.false.
1921         enddo
1922       enddo
1923 !
1924 ! Read side-chain interaction parameters.
1925 !
1926 !el from module energy - COMMON.INTERACT-------
1927       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1928       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1929       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1930       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1931       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1932       allocate(epslip(ntyp,ntyp))
1933       augm(:,:)=0.0D0
1934       chip(:)=0.0D0
1935       alp(:)=0.0D0
1936       sigma0(:)=0.0D0
1937       sigii(:)=0.0D0
1938       rr0(:)=0.0D0
1939
1940 !--------------------------------
1941
1942       read (isidep,*,end=117,err=117) ipot,expon
1943       if (ipot.lt.1 .or. ipot.gt.5) then
1944         write (iout,'(2a)') 'Error while reading SC interaction',&
1945                      'potential file - unknown potential type.'
1946 #ifdef MPI
1947         call MPI_Finalize(Ierror)
1948 #endif
1949         stop
1950       endif
1951       expon2=expon/2
1952       if(me.eq.king) &
1953        write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1954        ', exponents are ',expon,2*expon
1955 !      goto (10,20,30,30,40) ipot
1956       select case(ipot)
1957 !----------------------- LJ potential ---------------------------------
1958        case (1)
1959 !   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1960          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1961            (sigma0(i),i=1,ntyp)
1962         if (lprint) then
1963           write (iout,'(/a/)') 'Parameters of the LJ potential:'
1964           write (iout,'(a/)') 'The epsilon array:'
1965           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1966           write (iout,'(/a)') 'One-body parameters:'
1967           write (iout,'(a,4x,a)') 'residue','sigma'
1968           write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
1969         endif
1970 !      goto 50
1971 !----------------------- LJK potential --------------------------------
1972        case(2)
1973 !   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1974          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1975           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1976         if (lprint) then
1977           write (iout,'(/a/)') 'Parameters of the LJK potential:'
1978           write (iout,'(a/)') 'The epsilon array:'
1979           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1980           write (iout,'(/a)') 'One-body parameters:'
1981           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
1982           write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
1983                 i=1,ntyp)
1984         endif
1985 !      goto 50
1986 !---------------------- GB or BP potential -----------------------------
1987        case(3:4)
1988 !   30 do i=1,ntyp
1989         do i=1,ntyp
1990          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1991         enddo
1992         read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1993         read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1994         read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1995         read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1996         do i=1,ntyp
1997          read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
1998         enddo
1999
2000 ! For the GB potential convert sigma'**2 into chi'
2001         if (ipot.eq.4) then
2002           do i=1,ntyp
2003             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2004           enddo
2005         endif
2006         if (lprint) then
2007           write (iout,'(/a/)') 'Parameters of the BP potential:'
2008           write (iout,'(a/)') 'The epsilon array:'
2009           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2010           write (iout,'(/a)') 'One-body parameters:'
2011           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
2012                '    chip  ','    alph  '
2013           write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2014                              chip(i),alp(i),i=1,ntyp)
2015         endif
2016 !      goto 50
2017 !--------------------- GBV potential -----------------------------------
2018        case(5)
2019 !   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2020         read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2021           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2022           (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2023         if (lprint) then
2024           write (iout,'(/a/)') 'Parameters of the GBV potential:'
2025           write (iout,'(a/)') 'The epsilon array:'
2026           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2027           write (iout,'(/a)') 'One-body parameters:'
2028           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
2029               's||/s_|_^2','    chip  ','    alph  '
2030           write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2031                    sigii(i),chip(i),alp(i),i=1,ntyp)
2032         endif
2033        case default
2034         write(iout,*)"Wrong ipot"
2035 !   50 continue
2036       end select
2037       continue
2038       close (isidep)
2039 !-----------------------------------------------------------------------
2040 ! Calculate the "working" parameters of SC interactions.
2041
2042 !el from module energy - COMMON.INTERACT-------
2043       allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2044       allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2045       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2046       allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2047         dcavtub(ntyp1))
2048       allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2049         tubetranene(ntyp1))
2050       aa_aq(:,:)=0.0D0
2051       bb_aq(:,:)=0.0D0
2052       aa_lip(:,:)=0.0D0
2053       bb_lip(:,:)=0.0D0
2054       chi(:,:)=0.0D0
2055       sigma(:,:)=0.0D0
2056       r0(:,:)=0.0D0
2057       acavtub(:)=0.0d0
2058       bcavtub(:)=0.0d0
2059       ccavtub(:)=0.0d0
2060       dcavtub(:)=0.0d0
2061       sc_aa_tube_par(:)=0.0d0
2062       sc_bb_tube_par(:)=0.0d0
2063
2064 !--------------------------------
2065
2066       do i=2,ntyp
2067         do j=1,i-1
2068           eps(i,j)=eps(j,i)
2069           epslip(i,j)=epslip(j,i)
2070         enddo
2071       enddo
2072       do i=1,ntyp
2073         do j=i,ntyp
2074           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2075           sigma(j,i)=sigma(i,j)
2076           rs0(i,j)=dwa16*sigma(i,j)
2077           rs0(j,i)=rs0(i,j)
2078         enddo
2079       enddo
2080       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2081        'Working parameters of the SC interactions:',&
2082        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
2083        '  chi1   ','   chi2   '
2084       do i=1,ntyp
2085         do j=i,ntyp
2086           epsij=eps(i,j)
2087           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2088             rrij=sigma(i,j)
2089           else
2090             rrij=rr0(i)+rr0(j)
2091           endif
2092           r0(i,j)=rrij
2093           r0(j,i)=rrij
2094           rrij=rrij**expon
2095           epsij=eps(i,j)
2096           sigeps=dsign(1.0D0,epsij)
2097           epsij=dabs(epsij)
2098           aa_aq(i,j)=epsij*rrij*rrij
2099           bb_aq(i,j)=-sigeps*epsij*rrij
2100           aa_aq(j,i)=aa_aq(i,j)
2101           bb_aq(j,i)=bb_aq(i,j)
2102           epsijlip=epslip(i,j)
2103           sigeps=dsign(1.0D0,epsijlip)
2104           epsijlip=dabs(epsijlip)
2105           aa_lip(i,j)=epsijlip*rrij*rrij
2106           bb_lip(i,j)=-sigeps*epsijlip*rrij
2107           aa_lip(j,i)=aa_lip(i,j)
2108           bb_lip(j,i)=bb_lip(i,j)
2109 !C          write(iout,*) aa_lip
2110           if (ipot.gt.2) then
2111             sigt1sq=sigma0(i)**2
2112             sigt2sq=sigma0(j)**2
2113             sigii1=sigii(i)
2114             sigii2=sigii(j)
2115             ratsig1=sigt2sq/sigt1sq
2116             ratsig2=1.0D0/ratsig1
2117             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2118             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2119             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2120           else
2121             rsum_max=sigma(i,j)
2122           endif
2123 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2124             sigmaii(i,j)=rsum_max
2125             sigmaii(j,i)=rsum_max
2126 !         else
2127 !           sigmaii(i,j)=r0(i,j)
2128 !           sigmaii(j,i)=r0(i,j)
2129 !         endif
2130 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2131           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2132             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2133             augm(i,j)=epsij*r_augm**(2*expon)
2134 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2135             augm(j,i)=augm(i,j)
2136           else
2137             augm(i,j)=0.0D0
2138             augm(j,i)=0.0D0
2139           endif
2140           if (lprint) then
2141             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2142             restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2143             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2144           endif
2145         enddo
2146       enddo
2147       write(iout,*) "tube param"
2148       read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2149       ccavtubpep,dcavtubpep,tubetranenepep
2150       sigmapeptube=sigmapeptube**6
2151       sigeps=dsign(1.0D0,epspeptube)
2152       epspeptube=dabs(epspeptube)
2153       pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2154       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2155       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2156       do i=1,ntyp
2157        read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2158       ccavtub(i),dcavtub(i),tubetranene(i)
2159        sigmasctube=sigmasctube**6
2160        sigeps=dsign(1.0D0,epssctube)
2161        epssctube=dabs(epssctube)
2162        sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2163        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2164       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2165       enddo
2166
2167       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2168       bad(:,:)=0.0D0
2169
2170 #ifdef OLDSCP
2171 !
2172 ! Define the SC-p interaction constants (hard-coded; old style)
2173 !
2174       do i=1,ntyp
2175 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2176 ! helix formation)
2177 !       aad(i,1)=0.3D0*4.0D0**12
2178 ! Following line for constants currently implemented
2179 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2180         aad(i,1)=1.5D0*4.0D0**12
2181 !       aad(i,1)=0.17D0*5.6D0**12
2182         aad(i,2)=aad(i,1)
2183 ! "Soft" SC-p repulsion
2184         bad(i,1)=0.0D0
2185 ! Following line for constants currently implemented
2186 !       aad(i,1)=0.3D0*4.0D0**6
2187 ! "Hard" SC-p repulsion
2188         bad(i,1)=3.0D0*4.0D0**6
2189 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2190         bad(i,2)=bad(i,1)
2191 !       aad(i,1)=0.0D0
2192 !       aad(i,2)=0.0D0
2193 !       bad(i,1)=1228.8D0
2194 !       bad(i,2)=1228.8D0
2195       enddo
2196 #else
2197 !
2198 ! 8/9/01 Read the SC-p interaction constants from file
2199 !
2200       do i=1,ntyp
2201         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2202       enddo
2203       do i=1,ntyp
2204         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2205         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2206         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2207         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2208       enddo
2209 !      lprint=.true.
2210       if (lprint) then
2211         write (iout,*) "Parameters of SC-p interactions:"
2212         do i=1,ntyp
2213           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2214            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2215         enddo
2216       endif
2217 !      lprint=.false.
2218 #endif
2219 !
2220 ! Define the constants of the disulfide bridge
2221 !
2222       ebr=-5.50D0
2223 !
2224 ! Old arbitrary potential - commented out.
2225 !
2226 !      dbr= 4.20D0
2227 !      fbr= 3.30D0
2228 !
2229 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2230 ! energy surface of diethyl disulfide.
2231 ! A. Liwo and U. Kozlowska, 11/24/03
2232 !
2233       D0CM = 3.78d0
2234       AKCM = 15.1d0
2235       AKTH = 11.0d0
2236       AKCT = 12.0d0
2237       V1SS =-1.08d0
2238       V2SS = 7.61d0
2239       V3SS = 13.7d0
2240 !      akcm=0.0d0
2241 !      akth=0.0d0
2242 !      akct=0.0d0
2243 !      v1ss=0.0d0
2244 !      v2ss=0.0d0
2245 !      v3ss=0.0d0
2246
2247       if(me.eq.king) then
2248       write (iout,'(/a)') "Disulfide bridge parameters:"
2249       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2250       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2251       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2252       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2253         ' v3ss:',v3ss
2254       endif
2255       return
2256   111 write (iout,*) "Error reading bending energy parameters."
2257       goto 999
2258   112 write (iout,*) "Error reading rotamer energy parameters."
2259       goto 999
2260   113 write (iout,*) "Error reading torsional energy parameters."
2261       goto 999
2262   114 write (iout,*) "Error reading double torsional energy parameters."
2263       goto 999
2264   115 write (iout,*) &
2265         "Error reading cumulant (multibody energy) parameters."
2266       goto 999
2267   116 write (iout,*) "Error reading electrostatic energy parameters."
2268       goto 999
2269   117 write (iout,*) "Error reading side chain interaction parameters."
2270       goto 999
2271   118 write (iout,*) "Error reading SCp interaction parameters."
2272       goto 999
2273   119 write (iout,*) "Error reading SCCOR parameters"
2274   999 continue
2275 #ifdef MPI
2276       call MPI_Finalize(Ierror)
2277 #endif
2278       stop
2279       return
2280       end subroutine parmread
2281 #endif
2282 !-----------------------------------------------------------------------------
2283 ! printmat.f
2284 !-----------------------------------------------------------------------------
2285       subroutine printmat(ldim,m,n,iout,key,a)
2286
2287       integer :: n,ldim
2288       character(len=3),dimension(n) :: key
2289       real(kind=8),dimension(ldim,n) :: a
2290 !el local variables
2291       integer :: i,j,k,m,iout,nlim
2292
2293       do 1 i=1,n,8
2294       nlim=min0(i+7,n)
2295       write (iout,1000) (key(k),k=i,nlim)
2296       write (iout,1020)
2297  1000 format (/5x,8(6x,a3))
2298  1020 format (/80(1h-)/)
2299       do 2 j=1,n
2300       write (iout,1010) key(j),(a(j,k),k=i,nlim)
2301     2 continue
2302     1 continue
2303  1010 format (a3,2x,8(f9.4))
2304       return
2305       end subroutine printmat
2306 !-----------------------------------------------------------------------------
2307 ! readpdb.F
2308 !-----------------------------------------------------------------------------
2309       subroutine readpdb
2310 ! Read the PDB file and convert the peptide geometry into virtual-chain
2311 ! geometry.
2312       use geometry_data
2313       use energy_data, only: itype
2314       use control_data
2315       use compare_data
2316       use MPI_data
2317       use control, only: rescode,sugarcode
2318 !      implicit real*8 (a-h,o-z)
2319 !      include 'DIMENSIONS'
2320 !      include 'COMMON.LOCAL'
2321 !      include 'COMMON.VAR'
2322 !      include 'COMMON.CHAIN'
2323 !      include 'COMMON.INTERACT'
2324 !      include 'COMMON.IOUNITS'
2325 !      include 'COMMON.GEO'
2326 !      include 'COMMON.NAMES'
2327 !      include 'COMMON.CONTROL'
2328 !      include 'COMMON.DISTFIT'
2329 !      include 'COMMON.SETUP'
2330       integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2331 !        ishift_pdb
2332       logical :: lprn=.true.,fail
2333       real(kind=8),dimension(3) :: e1,e2,e3
2334       real(kind=8) :: dcj,efree_temp
2335       character(len=3) :: seq,res
2336       character(len=5) :: atom
2337       character(len=80) :: card
2338       real(kind=8),dimension(3,20) :: sccor
2339       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin       !rescode,
2340       integer :: isugar
2341       character*1 :: sugar
2342       real(kind=8) :: cou
2343       real(kind=8),dimension(3) :: ccc
2344 !el local varables
2345       integer,dimension(2,maxres/3) :: hfrag_alloc
2346       integer,dimension(4,maxres/3) :: bfrag_alloc
2347       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2348       real(kind=8),dimension(:,:), allocatable  :: c_temporary
2349       integer,dimension(:,:) , allocatable  :: itype_temporary
2350       efree_temp=0.0d0
2351       ibeg=1
2352       ishift1=0
2353       ishift=0
2354       molecule=1
2355       counter=0
2356 !      write (2,*) "UNRES_PDB",unres_pdb
2357       ires=0
2358       ires_old=0
2359       nres=0
2360       iii=0
2361       lsecondary=.false.
2362       nhfrag=0
2363       nbfrag=0
2364 !-----------------------------
2365       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2366       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2367
2368       do i=1,100000
2369         read (ipdbin,'(a80)',end=10) card
2370 !       write (iout,'(a)') card
2371         if (card(:5).eq.'HELIX') then
2372           nhfrag=nhfrag+1
2373           lsecondary=.true.
2374           read(card(22:25),*) hfrag(1,nhfrag)
2375           read(card(34:37),*) hfrag(2,nhfrag)
2376         endif
2377         if (card(:5).eq.'SHEET') then
2378           nbfrag=nbfrag+1
2379           lsecondary=.true.
2380           read(card(24:26),*) bfrag(1,nbfrag)
2381           read(card(35:37),*) bfrag(2,nbfrag)
2382 !rc----------------------------------------
2383 !rc  to be corrected !!!
2384           bfrag(3,nbfrag)=bfrag(1,nbfrag)
2385           bfrag(4,nbfrag)=bfrag(2,nbfrag)
2386 !rc----------------------------------------
2387         endif
2388         if (card(:3).eq.'END') then
2389           goto 10
2390         else if (card(:3).eq.'TER') then
2391 ! End current chain
2392           ires_old=ires+2
2393           ishift=ishift+1
2394           ishift1=ishift1+1
2395           itype(ires_old,molecule)=ntyp1_molec(molecule)
2396           itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2397           nres_molec(molecule)=nres_molec(molecule)+2
2398           ibeg=2
2399 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2400           if (unres_pdb) then
2401             do j=1,3
2402               dc(j,ires)=sccor(j,iii)
2403             enddo
2404           else
2405             call sccenter(ires,iii,sccor)
2406 !          iii=0
2407           endif
2408           iii=0
2409         endif
2410 ! Read free energy
2411         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2412 ! Fish out the ATOM cards.
2413 !        write(iout,*) 'card',card(1:20)
2414         if (index(card(1:4),'ATOM').gt.0) then
2415           read (card(12:16),*) atom
2416 !          write (iout,*) "! ",atom," !",ires
2417 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2418           read (card(23:26),*) ires
2419           read (card(18:20),'(a3)') res
2420 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2421 !     &      " ires_old",ires_old
2422 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2423 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2424           if (ires-ishift+ishift1.ne.ires_old) then
2425 ! Calculate the CM of the preceding residue.
2426 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2427             if (ibeg.eq.0) then
2428 !              write (iout,*) "Calculating sidechain center iii",iii
2429               if (unres_pdb) then
2430                 do j=1,3
2431                   dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2432                 enddo
2433               else
2434                 call sccenter(ires_old,iii,sccor)
2435               endif
2436               iii=0
2437             endif
2438 ! Start new residue.
2439             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2440               ires=ires_old
2441               cycle
2442             else if (ibeg.eq.1) then
2443               write (iout,*) "BEG ires",ires
2444               ishift=ires-1
2445               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2446                 ishift=ishift-1
2447                 itype(1,1)=ntyp1
2448                 nres_molec(molecule)=nres_molec(molecule)+1
2449               endif
2450               ires=ires-ishift+ishift1
2451               ires_old=ires
2452 !              write (iout,*) "ishift",ishift," ires",ires,&
2453 !               " ires_old",ires_old
2454               ibeg=0
2455             else if (ibeg.eq.2) then
2456 ! Start a new chain
2457               ishift=-ires_old+ires-1 !!!!!
2458               ishift1=ishift1-1    !!!!!
2459 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2460               ires=ires-ishift+ishift1
2461               print *,ires,ishift,ishift1
2462               ires_old=ires
2463               ibeg=0
2464             else
2465               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2466               ires=ires-ishift+ishift1
2467               ires_old=ires
2468             endif
2469 !            print *,'atom',ires,atom
2470             if (res.eq.'ACE' .or. res.eq.'NHE') then
2471               itype(ires,1)=10
2472             else
2473              if (atom.eq.'CA  '.or.atom.eq.'N   ') then
2474              molecule=1
2475               itype(ires,molecule)=rescode(ires,res,0,molecule)
2476 !              nres_molec(molecule)=nres_molec(molecule)+1
2477             else
2478              molecule=2
2479               itype(ires,molecule)=rescode(ires,res(3:4),0,molecule)
2480 !              nres_molec(molecule)=nres_molec(molecule)+1
2481             endif
2482             endif
2483           else
2484             ires=ires-ishift+ishift1
2485           endif
2486 !          write (iout,*) "ires_old",ires_old," ires",ires
2487           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2488 !            ishift1=ishift1+1
2489           endif
2490 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2491           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2492              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2493             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2494 !              print *,ires,ishift,ishift1
2495 !            write (iout,*) "backbone ",atom
2496 #ifdef DEBUG
2497             write (iout,'(2i3,2x,a,3f8.3)') &
2498             ires,itype(ires,1),res,(c(j,ires),j=1,3)
2499 #endif
2500             iii=iii+1
2501               nres_molec(molecule)=nres_molec(molecule)+1
2502             do j=1,3
2503               sccor(j,iii)=c(j,ires)
2504             enddo
2505           else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2506                atom.eq."C2'" .or. atom.eq."C3'" &
2507                .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2508             read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2509 !c            write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2510               print *,ires,ishift,ishift1
2511             counter=counter+1
2512 !            iii=iii+1
2513 !            do j=1,3
2514 !              sccor(j,iii)=c(j,ires)
2515 !            enddo
2516             do j=1,3
2517               c(j,ires)=c(j,ires)+ccc(j)/5.0
2518             enddo
2519              if (counter.eq.5) then
2520 !            iii=iii+1
2521               nres_molec(molecule)=nres_molec(molecule)+1
2522 !            do j=1,3
2523 !              sccor(j,iii)=c(j,ires)
2524 !            enddo
2525              counter=0
2526            endif
2527 !            print *, "ATOM",atom(1:3)
2528           else if (atom(1:3).eq."C5'") then
2529              read (card(19:19),'(a1)') sugar
2530              isugar=sugarcode(sugar,ires)
2531             if (ibeg.eq.1) then
2532               istype(1)=isugar
2533             else
2534               istype(ires)=isugar
2535             endif
2536             if (unres_pdb) then
2537               read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2538             else
2539               iii=iii+1
2540               read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2541             endif
2542 !            write (*,*) card(23:27),ires,itype(ires,1)
2543           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2544                    atom.ne.'N' .and. atom.ne.'C' .and. &
2545                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2546                    atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2547                    .and. atom.ne.'P  '.and. &
2548                   atom(1:1).ne.'H' .and. &
2549                   atom.ne.'OP1' .and. atom.ne.'OP2 ') then
2550 !            write (iout,*) "sidechain ",atom
2551 !            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2552                  if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
2553 !                        write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
2554             iii=iii+1
2555             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2556               endif
2557           endif
2558         endif
2559       enddo
2560    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2561       if (ires.eq.0) return
2562 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2563 ! system
2564       nres=ires
2565       if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
2566 !      print *,'I have', nres_molec(:)
2567
2568       do k=1,5
2569        if (nres_molec(k).eq.0) cycle
2570       do i=2,nres-1
2571 !        write (iout,*) i,itype(i,1)
2572 !        if (itype(i,1).eq.ntyp1) then
2573 !          write (iout,*) "dummy",i,itype(i,1)
2574 !          do j=1,3
2575 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2576 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
2577 !            dc(j,i)=c(j,i)
2578 !          enddo
2579 !        endif
2580         if (itype(i,k).eq.ntyp1_molec(k)) then
2581          if (itype(i+1,k).eq.ntyp1_molec(k)) then
2582           if (itype(i+2,k).eq.0) then
2583            print *,"masz sieczke"
2584            do j=1,5
2585             if (itype(i+2,j).ne.0) then
2586             itype(i+1,k)=0
2587             itype(i+1,j)=ntyp1_molec(j)
2588             nres_molec(k)=nres_molec(k)-1
2589             nres_molec(j)=nres_molec(j)+1
2590             go to 3331
2591             endif
2592            enddo
2593  3331      continue
2594           endif
2595 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2596 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
2597 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
2598            if (unres_pdb) then
2599 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2600 !            print *,i,'tu dochodze'
2601             call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2602             if (fail) then
2603               e2(1)=0.0d0
2604               e2(2)=1.0d0
2605               e2(3)=0.0d0
2606             endif !fail
2607             print *,i,'a tu?'
2608             do j=1,3
2609              c(j,i)=c(j,i-1)-1.9d0*e2(j)
2610             enddo
2611            else   !unres_pdb
2612            do j=1,3
2613              dcj=(c(j,i-2)-c(j,i-3))/2.0
2614             if (dcj.eq.0) dcj=1.23591524223
2615              c(j,i)=c(j,i-1)+dcj
2616              c(j,nres+i)=c(j,i)
2617            enddo
2618           endif   !unres_pdb
2619          else     !itype(i+1,1).eq.ntyp1
2620           if (unres_pdb) then
2621 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2622             call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2623             if (fail) then
2624               e2(1)=0.0d0
2625               e2(2)=1.0d0
2626               e2(3)=0.0d0
2627             endif
2628             do j=1,3
2629               c(j,i)=c(j,i+1)-1.9d0*e2(j)
2630             enddo
2631           else !unres_pdb
2632            do j=1,3
2633             dcj=(c(j,i+3)-c(j,i+2))/2.0
2634             if (dcj.eq.0) dcj=1.23591524223
2635             c(j,i)=c(j,i+1)-dcj
2636             c(j,nres+i)=c(j,i)
2637            enddo
2638           endif !unres_pdb
2639          endif !itype(i+1,1).eq.ntyp1
2640         endif  !itype.eq.ntyp1
2641
2642       enddo
2643       enddo
2644 ! Calculate the CM of the last side chain.
2645       if (iii.gt.0)  then
2646       if (unres_pdb) then
2647         do j=1,3
2648           dc(j,ires)=sccor(j,iii)
2649         enddo
2650       else
2651         call sccenter(ires,iii,sccor)
2652       endif
2653       endif
2654 !      nres=ires
2655       nsup=nres
2656       nstart_sup=1
2657 !      print *,"molecule",molecule
2658       if (itype(nres,1).ne.10) then
2659         nres=nres+1
2660         itype(nres,molecule)=ntyp1_molec(molecule)
2661         nres_molec(molecule)=nres_molec(molecule)+1
2662         if (unres_pdb) then
2663 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2664           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2665           if (fail) then
2666             e2(1)=0.0d0
2667             e2(2)=1.0d0
2668             e2(3)=0.0d0
2669           endif
2670           do j=1,3
2671             c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2672           enddo
2673         else
2674         do j=1,3
2675           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2676           c(j,nres)=c(j,nres-1)+dcj
2677           c(j,2*nres)=c(j,nres)
2678         enddo
2679         endif
2680       endif
2681 !     print *,'I have',nres, nres_molec(:)
2682
2683 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2684 #ifdef WHAM_RUN
2685       if (nres.ne.nres0) then
2686         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2687                        " NRES0=",nres0
2688         stop "Error nres value in WHAM input"
2689       endif
2690 #endif
2691 !---------------------------------
2692 !el reallocate tables
2693 !      do i=1,maxres/3
2694 !       do j=1,2
2695 !         hfrag_alloc(j,i)=hfrag(j,i)
2696 !        enddo
2697 !       do j=1,4
2698 !         bfrag_alloc(j,i)=bfrag(j,i)
2699 !        enddo
2700 !      enddo
2701
2702 !      deallocate(hfrag)
2703 !      deallocate(bfrag)
2704 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
2705 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2706 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2707 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
2708
2709 !      do i=1,nhfrag
2710 !       do j=1,2
2711 !         hfrag(j,i)=hfrag_alloc(j,i)
2712 !        enddo
2713 !      enddo
2714 !      do i=1,nbfrag
2715 !       do j=1,4
2716 !         bfrag(j,i)=bfrag_alloc(j,i)
2717 !        enddo
2718 !      enddo
2719 !el end reallocate tables
2720 !---------------------------------
2721       do i=2,nres-1
2722         do j=1,3
2723           c(j,i+nres)=dc(j,i)
2724         enddo
2725       enddo
2726       do j=1,3
2727         c(j,nres+1)=c(j,1)
2728         c(j,2*nres)=c(j,nres)
2729       enddo
2730
2731       if (itype(1,1).eq.ntyp1) then
2732         nsup=nsup-1
2733         nstart_sup=2
2734         if (unres_pdb) then
2735 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2736           call refsys(2,3,4,e1,e2,e3,fail)
2737           if (fail) then
2738             e2(1)=0.0d0
2739             e2(2)=1.0d0
2740             e2(3)=0.0d0
2741           endif
2742           do j=1,3
2743             c(j,1)=c(j,2)-1.9d0*e2(j)
2744           enddo
2745         else
2746         do j=1,3
2747           dcj=(c(j,4)-c(j,3))/2.0
2748           c(j,1)=c(j,2)-dcj
2749           c(j,nres+1)=c(j,1)
2750         enddo
2751         endif
2752       endif
2753 ! First lets assign correct dummy to correct type of chain
2754 ! 1) First residue
2755       if (itype(1,1).eq.ntyp1) then
2756         if (itype(2,1).eq.0) then
2757          do j=2,5
2758            if (itype(2,j).ne.0) then
2759            itype(1,1)=0
2760            itype(1,j)=ntyp1_molec(j)
2761            nres_molec(1)=nres_molec(1)-1
2762            nres_molec(j)=nres_molec(j)+1
2763            go to 3231
2764            endif
2765          enddo
2766 3231    continue
2767         endif
2768        endif
2769        print *,'I have',nres, nres_molec(:)
2770
2771 ! Copy the coordinates to reference coordinates
2772 !      do i=1,2*nres
2773 !        do j=1,3
2774 !          cref(j,i)=c(j,i)
2775 !        enddo
2776 !      enddo
2777 ! Calculate internal coordinates.
2778       if (lprn) then
2779       write (iout,'(/a)') &
2780         "Cartesian coordinates of the reference structure"
2781       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2782        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2783       do ires=1,nres
2784         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
2785           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
2786           (c(j,ires+nres),j=1,3)
2787       enddo
2788       endif
2789 ! znamy już nres wiec mozna alokowac tablice
2790 ! Calculate internal coordinates.
2791       if(me.eq.king.or..not.out1file)then
2792        write (iout,'(a)') &
2793          "Backbone and SC coordinates as read from the PDB"
2794        do ires=1,nres
2795         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2796           ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
2797           (c(j,nres+ires),j=1,3)
2798        enddo
2799       endif
2800 ! NOW LETS ROCK! SORTING
2801       allocate(c_temporary(3,2*nres))
2802       allocate(itype_temporary(nres,5))
2803        itype_temporary(:,:)=0
2804       seqalingbegin=1
2805       do k=1,5
2806         do i=1,nres
2807          if (itype(i,k).ne.0) then
2808           do j=1,3
2809           c_temporary(j,seqalingbegin)=c(j,i)
2810           c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
2811
2812           enddo
2813           itype_temporary(seqalingbegin,k)=itype(i,k)
2814           seqalingbegin=seqalingbegin+1
2815          endif
2816         enddo
2817        enddo
2818        do i=1,2*nres
2819         do j=1,3
2820         c(j,i)=c_temporary(j,i)
2821         enddo
2822        enddo
2823        do k=1,5
2824         do i=1,nres
2825          itype(i,k)=itype_temporary(i,k)
2826         enddo
2827        enddo
2828       if (itype(1,1).eq.ntyp1) then
2829         nsup=nsup-1
2830         nstart_sup=2
2831         if (unres_pdb) then
2832 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2833           call refsys(2,3,4,e1,e2,e3,fail)
2834           if (fail) then
2835             e2(1)=0.0d0
2836             e2(2)=1.0d0
2837             e2(3)=0.0d0
2838           endif
2839           do j=1,3
2840             c(j,1)=c(j,2)-1.9d0*e2(j)
2841           enddo
2842         else
2843         do j=1,3
2844           dcj=(c(j,4)-c(j,3))/2.0
2845           c(j,1)=c(j,2)-dcj
2846           c(j,nres+1)=c(j,1)
2847         enddo
2848         endif
2849       endif
2850
2851       if (lprn) then
2852       write (iout,'(/a)') &
2853         "Cartesian coordinates of the reference structure after sorting"
2854       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2855        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2856       do ires=1,nres
2857         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
2858           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
2859           (c(j,ires+nres),j=1,3)
2860       enddo
2861       endif
2862
2863        print *,seqalingbegin,nres
2864       if(.not.allocated(vbld)) then
2865        allocate(vbld(2*nres))
2866        do i=1,2*nres
2867          vbld(i)=0.d0
2868        enddo
2869       endif
2870       if(.not.allocated(vbld_inv)) then
2871        allocate(vbld_inv(2*nres))
2872        do i=1,2*nres
2873          vbld_inv(i)=0.d0
2874        enddo
2875       endif
2876 !!!el
2877       if(.not.allocated(theta)) then
2878         allocate(theta(nres+2))
2879         theta(:)=0.0d0
2880       endif
2881
2882       if(.not.allocated(phi)) allocate(phi(nres+2))
2883       if(.not.allocated(alph)) allocate(alph(nres+2))
2884       if(.not.allocated(omeg)) allocate(omeg(nres+2))
2885       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2886       if(.not.allocated(phiref)) allocate(phiref(nres+2))
2887       if(.not.allocated(costtab)) allocate(costtab(nres))
2888       if(.not.allocated(sinttab)) allocate(sinttab(nres))
2889       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2890       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2891       if(.not.allocated(xxref)) allocate(xxref(nres))
2892       if(.not.allocated(yyref)) allocate(yyref(nres))
2893       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2894       if(.not.allocated(dc_norm)) then
2895 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2896         allocate(dc_norm(3,0:2*nres+2))
2897         dc_norm(:,:)=0.d0
2898       endif
2899
2900       call int_from_cart(.true.,.false.)
2901       call sc_loc_geom(.false.)
2902       do i=1,nres
2903         thetaref(i)=theta(i)
2904         phiref(i)=phi(i)
2905       enddo
2906 !      do i=1,2*nres
2907 !        vbld_inv(i)=0.d0
2908 !        vbld(i)=0.d0
2909 !      enddo
2910
2911       do i=1,nres-1
2912         do j=1,3
2913           dc(j,i)=c(j,i+1)-c(j,i)
2914           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2915         enddo
2916       enddo
2917       do i=2,nres-1
2918         do j=1,3
2919           dc(j,i+nres)=c(j,i+nres)-c(j,i)
2920           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2921         enddo
2922 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2923 !        vbld_inv(i+nres)
2924       enddo
2925 !      call chainbuild
2926 ! Copy the coordinates to reference coordinates
2927 ! Splits to single chain if occurs
2928
2929 !      do i=1,2*nres
2930 !        do j=1,3
2931 !          cref(j,i,cou)=c(j,i)
2932 !        enddo
2933 !      enddo
2934 !
2935       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2936       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2937       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2938 !-----------------------------
2939       kkk=1
2940       lll=0
2941       cou=1
2942       do i=1,nres
2943       lll=lll+1
2944 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2945       if (i.gt.1) then
2946       if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
2947       chain_length=lll-1
2948       kkk=kkk+1
2949 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2950       lll=1
2951       endif
2952       endif
2953         do j=1,3
2954           cref(j,i,cou)=c(j,i)
2955           cref(j,i+nres,cou)=c(j,i+nres)
2956           if (i.le.nres) then
2957           chain_rep(j,lll,kkk)=c(j,i)
2958           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2959           endif
2960          enddo
2961       enddo
2962       write (iout,*) chain_length
2963       if (chain_length.eq.0) chain_length=nres
2964       do j=1,3
2965       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2966       chain_rep(j,chain_length+nres,symetr) &
2967       =chain_rep(j,chain_length+nres,1)
2968       enddo
2969 ! diagnostic
2970 !       write (iout,*) "spraw lancuchy",chain_length,symetr
2971 !       do i=1,4
2972 !         do kkk=1,chain_length
2973 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2974 !         enddo
2975 !        enddo
2976 ! enddiagnostic
2977 ! makes copy of chains
2978         write (iout,*) "symetr", symetr
2979       do j=1,3
2980       dc(j,0)=c(j,1)
2981       enddo
2982
2983       if (symetr.gt.1) then
2984        call permut(symetr)
2985        nperm=1
2986        do i=1,symetr
2987        nperm=nperm*i
2988        enddo
2989        do i=1,nperm
2990        write(iout,*) (tabperm(i,kkk),kkk=1,4)
2991        enddo
2992        do i=1,nperm
2993         cou=0
2994         do kkk=1,symetr
2995          icha=tabperm(i,kkk)
2996 !         write (iout,*) i,icha
2997          do lll=1,chain_length
2998           cou=cou+1
2999            if (cou.le.nres) then
3000            do j=1,3
3001             kupa=mod(lll,chain_length)
3002             iprzes=(kkk-1)*chain_length+lll
3003             if (kupa.eq.0) kupa=chain_length
3004 !            write (iout,*) "kupa", kupa
3005             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3006             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3007           enddo
3008           endif
3009          enddo
3010         enddo
3011        enddo
3012        endif
3013 !-koniec robienia kopii
3014 ! diag
3015       do kkk=1,nperm
3016       write (iout,*) "nowa struktura", nperm
3017       do i=1,nres
3018         write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3019       cref(2,i,kkk),&
3020       cref(3,i,kkk),cref(1,nres+i,kkk),&
3021       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3022       enddo
3023   100 format (//'              alpha-carbon coordinates       ',&
3024                 '     centroid coordinates'/ &
3025                 '       ', 6X,'X',11X,'Y',11X,'Z', &
3026                                 10X,'X',11X,'Y',11X,'Z')
3027   110 format (a,'(',i3,')',6f12.5)
3028
3029       enddo
3030 !c enddiag
3031       do j=1,nbfrag
3032         do i=1,4
3033          bfrag(i,j)=bfrag(i,j)-ishift
3034         enddo
3035       enddo
3036
3037       do j=1,nhfrag
3038         do i=1,2
3039          hfrag(i,j)=hfrag(i,j)-ishift
3040         enddo
3041       enddo
3042       ishift_pdb=ishift
3043
3044       return
3045       end subroutine readpdb
3046 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3047 !-----------------------------------------------------------------------------
3048 ! readrtns_CSA.F
3049 !-----------------------------------------------------------------------------
3050       subroutine read_control
3051 !
3052 ! Read contorl data
3053 !
3054 !      use geometry_data
3055       use comm_machsw
3056       use energy_data
3057       use control_data
3058       use compare_data
3059       use MCM_data
3060       use map_data
3061       use csa_data
3062       use MD_data
3063       use MPI_data
3064       use random, only: random_init
3065 !      implicit real*8 (a-h,o-z)
3066 !      include 'DIMENSIONS'
3067 #ifdef MP
3068       use prng, only:prng_restart
3069       include 'mpif.h'
3070       logical :: OKRandom!, prng_restart
3071       real(kind=8) :: r1
3072 #endif
3073 !      include 'COMMON.IOUNITS'
3074 !      include 'COMMON.TIME1'
3075 !      include 'COMMON.THREAD'
3076 !      include 'COMMON.SBRIDGE'
3077 !      include 'COMMON.CONTROL'
3078 !      include 'COMMON.MCM'
3079 !      include 'COMMON.MAP'
3080 !      include 'COMMON.HEADER'
3081 !      include 'COMMON.CSA'
3082 !      include 'COMMON.CHAIN'
3083 !      include 'COMMON.MUCA'
3084 !      include 'COMMON.MD'
3085 !      include 'COMMON.FFIELD'
3086 !      include 'COMMON.INTERACT'
3087 !      include 'COMMON.SETUP'
3088 !el      integer :: KDIAG,ICORFL,IXDR
3089 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3090       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3091         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
3092 !      character(len=80) :: ucase
3093       character(len=640) :: controlcard
3094
3095       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3096       integer i
3097
3098       nglob_csa=0
3099       eglob_csa=1d99
3100       nmin_csa=0
3101       read (INP,'(a)') titel
3102       call card_concat(controlcard,.true.)
3103 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3104 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3105       call reada(controlcard,'SEED',seed,0.0D0)
3106       call random_init(seed)
3107 ! Set up the time limit (caution! The time must be input in minutes!)
3108       read_cart=index(controlcard,'READ_CART').gt.0
3109       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3110       call readi(controlcard,'SYM',symetr,1)
3111       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3112       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3113       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3114       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3115       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3116       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3117       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3118       call reada(controlcard,'DRMS',drms,0.1D0)
3119       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3120        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3121        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3122        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3123        write (iout,'(a,f10.1)')'DRMS    = ',drms
3124        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3125        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3126       endif
3127       call readi(controlcard,'NZ_START',nz_start,0)
3128       call readi(controlcard,'NZ_END',nz_end,0)
3129       call readi(controlcard,'IZ_SC',iz_sc,0)
3130       timlim=60.0D0*timlim
3131       safety = 60.0d0*safety
3132       timem=timlim
3133       modecalc=0
3134       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3135 !C SHIELD keyword sets if the shielding effect of side-chains is used
3136 !C 0 denots no shielding is used all peptide are equally despite the
3137 !C solvent accesible area
3138 !C 1 the newly introduced function
3139 !C 2 reseved for further possible developement
3140       call readi(controlcard,'SHIELD',shield_mode,0)
3141 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3142         write(iout,*) "shield_mode",shield_mode
3143 !C  Varibles set size of box
3144       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3145       protein=index(controlcard,"PROTEIN").gt.0
3146       ions=index(controlcard,"IONS").gt.0
3147       nucleic=index(controlcard,"NUCLEIC").gt.0
3148       write (iout,*) "with_theta_constr ",with_theta_constr
3149       AFMlog=(index(controlcard,'AFM'))
3150       selfguide=(index(controlcard,'SELFGUIDE'))
3151       print *,'AFMlog',AFMlog,selfguide,"KUPA"
3152       call readi(controlcard,'GENCONSTR',genconstr,0)
3153       call reada(controlcard,'BOXX',boxxsize,100.0d0)
3154       call reada(controlcard,'BOXY',boxysize,100.0d0)
3155       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3156       call readi(controlcard,'TUBEMOD',tubemode,0)
3157       write (iout,*) TUBEmode,"TUBEMODE"
3158       if (TUBEmode.gt.0) then
3159        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3160        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3161        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3162        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3163        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3164        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3165        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3166        buftubebot=bordtubebot+tubebufthick
3167        buftubetop=bordtubetop-tubebufthick
3168       endif
3169
3170 ! CUTOFFF ON ELECTROSTATICS
3171       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3172       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3173       write(iout,*) "R_CUT_ELE=",r_cut_ele
3174 ! Lipidic parameters
3175       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3176       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3177       if (lipthick.gt.0.0d0) then
3178        bordliptop=(boxzsize+lipthick)/2.0
3179        bordlipbot=bordliptop-lipthick
3180       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3181       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3182       buflipbot=bordlipbot+lipbufthick
3183       bufliptop=bordliptop-lipbufthick
3184       if ((lipbufthick*2.0d0).gt.lipthick) &
3185        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3186       endif !lipthick.gt.0
3187       write(iout,*) "bordliptop=",bordliptop
3188       write(iout,*) "bordlipbot=",bordlipbot
3189       write(iout,*) "bufliptop=",bufliptop
3190       write(iout,*) "buflipbot=",buflipbot
3191       write (iout,*) "SHIELD MODE",shield_mode
3192
3193 !C-------------------------
3194       minim=(index(controlcard,'MINIMIZE').gt.0)
3195       dccart=(index(controlcard,'CART').gt.0)
3196       overlapsc=(index(controlcard,'OVERLAP').gt.0)
3197       overlapsc=.not.overlapsc
3198       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3199       searchsc=.not.searchsc
3200       sideadd=(index(controlcard,'SIDEADD').gt.0)
3201       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3202       outpdb=(index(controlcard,'PDBOUT').gt.0)
3203       outmol2=(index(controlcard,'MOL2OUT').gt.0)
3204       pdbref=(index(controlcard,'PDBREF').gt.0)
3205       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3206       indpdb=index(controlcard,'PDBSTART')
3207       extconf=(index(controlcard,'EXTCONF').gt.0)
3208       call readi(controlcard,'IPRINT',iprint,0)
3209       call readi(controlcard,'MAXGEN',maxgen,10000)
3210       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3211       call readi(controlcard,"KDIAG",kdiag,0)
3212       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3213       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3214        write (iout,*) "RESCALE_MODE",rescale_mode
3215       split_ene=index(controlcard,'SPLIT_ENE').gt.0
3216       if (index(controlcard,'REGULAR').gt.0.0D0) then
3217         call reada(controlcard,'WEIDIS',weidis,0.1D0)
3218         modecalc=1
3219         refstr=.true.
3220       endif
3221       if (index(controlcard,'CHECKGRAD').gt.0) then
3222         modecalc=5
3223         if (index(controlcard,'CART').gt.0) then
3224           icheckgrad=1
3225         elseif (index(controlcard,'CARINT').gt.0) then
3226           icheckgrad=2
3227         else
3228           icheckgrad=3
3229         endif
3230       elseif (index(controlcard,'THREAD').gt.0) then
3231         modecalc=2
3232         call readi(controlcard,'THREAD',nthread,0)
3233         if (nthread.gt.0) then
3234           call reada(controlcard,'WEIDIS',weidis,0.1D0)
3235         else
3236           if (fg_rank.eq.0) &
3237           write (iout,'(a)')'A number has to follow the THREAD keyword.'
3238           stop 'Error termination in Read_Control.'
3239         endif
3240       else if (index(controlcard,'MCMA').gt.0) then
3241         modecalc=3
3242       else if (index(controlcard,'MCEE').gt.0) then
3243         modecalc=6
3244       else if (index(controlcard,'MULTCONF').gt.0) then
3245         modecalc=4
3246       else if (index(controlcard,'MAP').gt.0) then
3247         modecalc=7
3248         call readi(controlcard,'MAP',nmap,0)
3249       else if (index(controlcard,'CSA').gt.0) then
3250         modecalc=8
3251 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3252 !rc
3253 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3254 !rc
3255 !rc        modecalc=9
3256 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3257 !fmc        modecalc=10
3258       else if (index(controlcard,'SOFTREG').gt.0) then
3259         modecalc=11
3260       else if (index(controlcard,'CHECK_BOND').gt.0) then
3261         modecalc=-1
3262       else if (index(controlcard,'TEST').gt.0) then
3263         modecalc=-2
3264       else if (index(controlcard,'MD').gt.0) then
3265         modecalc=12
3266       else if (index(controlcard,'RE ').gt.0) then
3267         modecalc=14
3268       endif
3269
3270       lmuca=index(controlcard,'MUCA').gt.0
3271       call readi(controlcard,'MUCADYN',mucadyn,0)
3272       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3273       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3274        then
3275        write (iout,*) 'MUCADYN=',mucadyn
3276        write (iout,*) 'MUCASMOOTH=',muca_smooth
3277       endif
3278
3279       iscode=index(controlcard,'ONE_LETTER')
3280       indphi=index(controlcard,'PHI')
3281       indback=index(controlcard,'BACK')
3282       iranconf=index(controlcard,'RAND_CONF')
3283       i2ndstr=index(controlcard,'USE_SEC_PRED')
3284       gradout=index(controlcard,'GRADOUT').gt.0
3285       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3286       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3287       if (me.eq.king .or. .not.out1file ) &
3288         write (iout,*) "DISTCHAINMAX",distchainmax
3289
3290       if(me.eq.king.or..not.out1file) &
3291        write (iout,'(2a)') diagmeth(kdiag),&
3292         ' routine used to diagonalize matrices.'
3293       if (shield_mode.gt.0) then
3294       pi=3.141592d0
3295 !C VSolvSphere the volume of solving sphere
3296 !C      print *,pi,"pi"
3297 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3298 !C there will be no distinction between proline peptide group and normal peptide
3299 !C group in case of shielding parameters
3300       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3301       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3302       write (iout,*) VSolvSphere,VSolvSphere_div
3303 !C long axis of side chain
3304       do i=1,ntyp
3305       long_r_sidechain(i)=vbldsc0(1,i)
3306       short_r_sidechain(i)=sigma0(i)
3307       write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3308          sigma0(i)
3309       enddo
3310       buff_shield=1.0d0
3311       endif
3312       return
3313       end subroutine read_control
3314 !-----------------------------------------------------------------------------
3315       subroutine read_REMDpar
3316 !
3317 ! Read REMD settings
3318 !
3319 !       use control
3320 !       use energy
3321 !       use geometry
3322       use REMD_data
3323       use MPI_data
3324       use control_data, only:out1file
3325 !      implicit real*8 (a-h,o-z)
3326 !      include 'DIMENSIONS'
3327 !      include 'COMMON.IOUNITS'
3328 !      include 'COMMON.TIME1'
3329 !      include 'COMMON.MD'
3330       use MD_data
3331 !el #ifndef LANG0
3332 !el      include 'COMMON.LANGEVIN'
3333 !el #else
3334 !el      include 'COMMON.LANGEVIN.lang0'
3335 !el #endif
3336 !      include 'COMMON.INTERACT'
3337 !      include 'COMMON.NAMES'
3338 !      include 'COMMON.GEO'
3339 !      include 'COMMON.REMD'
3340 !      include 'COMMON.CONTROL'
3341 !      include 'COMMON.SETUP'
3342 !      character(len=80) :: ucase
3343       character(len=320) :: controlcard
3344       character(len=3200) :: controlcard1
3345       integer :: iremd_m_total
3346 !el local variables
3347       integer :: i
3348 !     real(kind=8) :: var,ene
3349
3350       if(me.eq.king.or..not.out1file) &
3351        write (iout,*) "REMD setup"
3352
3353       call card_concat(controlcard,.true.)
3354       call readi(controlcard,"NREP",nrep,3)
3355       call readi(controlcard,"NSTEX",nstex,1000)
3356       call reada(controlcard,"RETMIN",retmin,10.0d0)
3357       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3358       mremdsync=(index(controlcard,'SYNC').gt.0)
3359       call readi(controlcard,"NSYN",i_sync_step,100)
3360       restart1file=(index(controlcard,'REST1FILE').gt.0)
3361       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3362       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3363       if(max_cache_traj_use.gt.max_cache_traj) &
3364                 max_cache_traj_use=max_cache_traj
3365       if(me.eq.king.or..not.out1file) then
3366 !d       if (traj1file) then
3367 !rc caching is in testing - NTWX is not ignored
3368 !d        write (iout,*) "NTWX value is ignored"
3369 !d        write (iout,*) "  trajectory is stored to one file by master"
3370 !d        write (iout,*) "  before exchange at NSTEX intervals"
3371 !d       endif
3372        write (iout,*) "NREP= ",nrep
3373        write (iout,*) "NSTEX= ",nstex
3374        write (iout,*) "SYNC= ",mremdsync
3375        write (iout,*) "NSYN= ",i_sync_step
3376        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3377       endif
3378       remd_tlist=.false.
3379       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3380       if (index(controlcard,'TLIST').gt.0) then
3381          remd_tlist=.true.
3382          call card_concat(controlcard1,.true.)
3383          read(controlcard1,*) (remd_t(i),i=1,nrep)
3384          if(me.eq.king.or..not.out1file) &
3385           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3386       endif
3387       remd_mlist=.false.
3388       if (index(controlcard,'MLIST').gt.0) then
3389          remd_mlist=.true.
3390          call card_concat(controlcard1,.true.)
3391          read(controlcard1,*) (remd_m(i),i=1,nrep)
3392          if(me.eq.king.or..not.out1file) then
3393           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3394           iremd_m_total=0
3395           do i=1,nrep
3396            iremd_m_total=iremd_m_total+remd_m(i)
3397           enddo
3398           write (iout,*) 'Total number of replicas ',iremd_m_total
3399          endif
3400       endif
3401       if(me.eq.king.or..not.out1file) &
3402        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3403       return
3404       end subroutine read_REMDpar
3405 !-----------------------------------------------------------------------------
3406       subroutine read_MDpar
3407 !
3408 ! Read MD settings
3409 !
3410       use control_data, only: r_cut,rlamb,out1file
3411       use energy_data
3412       use geometry_data, only: pi
3413       use MPI_data
3414 !      implicit real*8 (a-h,o-z)
3415 !      include 'DIMENSIONS'
3416 !      include 'COMMON.IOUNITS'
3417 !      include 'COMMON.TIME1'
3418 !      include 'COMMON.MD'
3419       use MD_data
3420 !el #ifndef LANG0
3421 !el      include 'COMMON.LANGEVIN'
3422 !el #else
3423 !el      include 'COMMON.LANGEVIN.lang0'
3424 !el #endif
3425 !      include 'COMMON.INTERACT'
3426 !      include 'COMMON.NAMES'
3427 !      include 'COMMON.GEO'
3428 !      include 'COMMON.SETUP'
3429 !      include 'COMMON.CONTROL'
3430 !      include 'COMMON.SPLITELE'
3431 !      character(len=80) :: ucase
3432       character(len=320) :: controlcard
3433 !el local variables
3434       integer :: i
3435       real(kind=8) :: eta
3436
3437       call card_concat(controlcard,.true.)
3438       call readi(controlcard,"NSTEP",n_timestep,1000000)
3439       call readi(controlcard,"NTWE",ntwe,100)
3440       call readi(controlcard,"NTWX",ntwx,1000)
3441       call reada(controlcard,"DT",d_time,1.0d-1)
3442       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3443       call reada(controlcard,"DAMAX",damax,1.0d1)
3444       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3445       call readi(controlcard,"LANG",lang,0)
3446       RESPA = index(controlcard,"RESPA") .gt. 0
3447       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3448       ntime_split0=ntime_split
3449       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3450       ntime_split0=ntime_split
3451       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3452       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3453       rest = index(controlcard,"REST").gt.0
3454       tbf = index(controlcard,"TBF").gt.0
3455       usampl = index(controlcard,"USAMPL").gt.0
3456       mdpdb = index(controlcard,"MDPDB").gt.0
3457       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3458       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3459       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3460       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3461       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3462       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3463       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3464       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3465       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3466       large = index(controlcard,"LARGE").gt.0
3467       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3468       rattle = index(controlcard,"RATTLE").gt.0
3469 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3470       nset=0
3471       if(usampl) then
3472         call read_fragments
3473       endif
3474
3475       if(me.eq.king.or..not.out1file) then
3476        write (iout,*)
3477        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3478        write (iout,*)
3479        write (iout,'(a)') "The units are:"
3480        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3481        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3482         " acceleration: angstrom/(48.9 fs)**2"
3483        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3484        write (iout,*)
3485        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3486        write (iout,'(a60,f10.5,a)') &
3487         "Initial time step of numerical integration:",d_time,&
3488         " natural units"
3489        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3490        if (RESPA) then
3491         write (iout,'(2a,i4,a)') &
3492           "A-MTS algorithm used; initial time step for fast-varying",&
3493           " short-range forces split into",ntime_split," steps."
3494         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3495          r_cut," lambda",rlamb
3496        endif
3497        write (iout,'(2a,f10.5)') &
3498         "Maximum acceleration threshold to reduce the time step",&
3499         "/increase split number:",damax
3500        write (iout,'(2a,f10.5)') &
3501         "Maximum predicted energy drift to reduce the timestep",&
3502         "/increase split number:",edriftmax
3503        write (iout,'(a60,f10.5)') &
3504        "Maximum velocity threshold to reduce velocities:",dvmax
3505        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3506        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3507        if (rattle) write (iout,'(a60)') &
3508         "Rattle algorithm used to constrain the virtual bonds"
3509       endif
3510       reset_fricmat=1000
3511       if (lang.gt.0) then
3512         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3513         call reada(controlcard,"RWAT",rwat,1.4d0)
3514         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3515         surfarea=index(controlcard,"SURFAREA").gt.0
3516         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3517         if(me.eq.king.or..not.out1file)then
3518          write (iout,'(/a,$)') "Langevin dynamics calculation"
3519          if (lang.eq.1) then
3520           write (iout,'(a/)') &
3521             " with direct integration of Langevin equations"
3522          else if (lang.eq.2) then
3523           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3524          else if (lang.eq.3) then
3525           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3526          else if (lang.eq.4) then
3527           write (iout,'(a/)') " in overdamped mode"
3528          else
3529           write (iout,'(//a,i5)') &
3530             "=========== ERROR: Unknown Langevin dynamics mode:",lang
3531           stop
3532          endif
3533          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3534          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3535          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3536          write (iout,'(a60,f10.5)') &
3537          "Scaling factor of the friction forces:",scal_fric
3538          if (surfarea) write (iout,'(2a,i10,a)') &
3539            "Friction coefficients will be scaled by solvent-accessible",&
3540            " surface area every",reset_fricmat," steps."
3541         endif
3542 ! Calculate friction coefficients and bounds of stochastic forces
3543         eta=6*pi*cPoise*etawat
3544         if(me.eq.king.or..not.out1file) &
3545          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3546           eta
3547         gamp=scal_fric*(pstok+rwat)*eta
3548         stdfp=dsqrt(2*Rb*t_bath/d_time)
3549         allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3550         do i=1,ntyp
3551           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3552           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3553         enddo
3554         if(me.eq.king.or..not.out1file)then
3555          write (iout,'(/2a/)') &
3556          "Radii of site types and friction coefficients and std's of",&
3557          " stochastic forces of fully exposed sites"
3558          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3559          do i=1,ntyp
3560           write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i),&
3561            gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3562          enddo
3563         endif
3564       else if (tbf) then
3565         if(me.eq.king.or..not.out1file)then
3566          write (iout,'(a)') "Berendsen bath calculation"
3567          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3568          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3569          if (reset_moment) &
3570          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3571          count_reset_moment," steps"
3572          if (reset_vel) &
3573           write (iout,'(a,i10,a)') &
3574           "Velocities will be reset at random every",count_reset_vel,&
3575          " steps"
3576         endif
3577       else
3578         if(me.eq.king.or..not.out1file) &
3579          write (iout,'(a31)') "Microcanonical mode calculation"
3580       endif
3581       if(me.eq.king.or..not.out1file)then
3582        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3583        if (usampl) then
3584           write(iout,*) "MD running with constraints."
3585           write(iout,*) "Equilibration time ", eq_time, " mtus."
3586           write(iout,*) "Constraining ", nfrag," fragments."
3587           write(iout,*) "Length of each fragment, weight and q0:"
3588           do iset=1,nset
3589            write (iout,*) "Set of restraints #",iset
3590            do i=1,nfrag
3591               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3592                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3593            enddo
3594            write(iout,*) "constraints between ", npair, "fragments."
3595            write(iout,*) "constraint pairs, weights and q0:"
3596            do i=1,npair
3597             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3598                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3599            enddo
3600            write(iout,*) "angle constraints within ", nfrag_back,&
3601             "backbone fragments."
3602            write(iout,*) "fragment, weights:"
3603            do i=1,nfrag_back
3604             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3605                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3606                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3607            enddo
3608           enddo
3609         iset=mod(kolor,nset)+1
3610        endif
3611       endif
3612       if(me.eq.king.or..not.out1file) &
3613        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3614       return
3615       end subroutine read_MDpar
3616 !-----------------------------------------------------------------------------
3617       subroutine map_read
3618
3619       use map_data
3620 !      implicit real*8 (a-h,o-z)
3621 !      include 'DIMENSIONS'
3622 !      include 'COMMON.MAP'
3623 !      include 'COMMON.IOUNITS'
3624       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3625       character(len=80) :: mapcard      !,ucase
3626 !el local variables
3627       integer :: imap
3628 !     real(kind=8) :: var,ene
3629
3630       do imap=1,nmap
3631         read (inp,'(a)') mapcard
3632         mapcard=ucase(mapcard)
3633         if (index(mapcard,'PHI').gt.0) then
3634           kang(imap)=1
3635         else if (index(mapcard,'THE').gt.0) then
3636           kang(imap)=2
3637         else if (index(mapcard,'ALP').gt.0) then
3638           kang(imap)=3
3639         else if (index(mapcard,'OME').gt.0) then
3640           kang(imap)=4
3641         else
3642           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3643           stop 'Error - illegal variable spec in MAP card.'
3644         endif
3645         call readi (mapcard,'RES1',res1(imap),0)
3646         call readi (mapcard,'RES2',res2(imap),0)
3647         if (res1(imap).eq.0) then
3648           res1(imap)=res2(imap)
3649         else if (res2(imap).eq.0) then
3650           res2(imap)=res1(imap)
3651         endif
3652         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3653           write (iout,'(a)') &
3654           'Error - illegal definition of variable group in MAP.'
3655           stop 'Error - illegal definition of variable group in MAP.'
3656         endif
3657         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3658         call reada(mapcard,'TO',ang_to(imap),0.0D0)
3659         call readi(mapcard,'NSTEP',nstep(imap),0)
3660         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3661           write (iout,'(a)') &
3662            'Illegal boundary and/or step size specification in MAP.'
3663           stop 'Illegal boundary and/or step size specification in MAP.'
3664         endif
3665       enddo ! imap
3666       return
3667       end subroutine map_read
3668 !-----------------------------------------------------------------------------
3669       subroutine csaread
3670
3671       use control_data, only: vdisulf
3672       use csa_data
3673 !      implicit real*8 (a-h,o-z)
3674 !      include 'DIMENSIONS'
3675 !      include 'COMMON.IOUNITS'
3676 !      include 'COMMON.GEO'
3677 !      include 'COMMON.CSA'
3678 !      include 'COMMON.BANK'
3679 !      include 'COMMON.CONTROL'
3680 !      character(len=80) :: ucase
3681       character(len=620) :: mcmcard
3682 !el local variables
3683 !     integer :: ntf,ik,iw_pdb
3684 !     real(kind=8) :: var,ene
3685
3686       call card_concat(mcmcard,.true.)
3687
3688       call readi(mcmcard,'NCONF',nconf,50)
3689       call readi(mcmcard,'NADD',nadd,0)
3690       call readi(mcmcard,'JSTART',jstart,1)
3691       call readi(mcmcard,'JEND',jend,1)
3692       call readi(mcmcard,'NSTMAX',nstmax,500000)
3693       call readi(mcmcard,'N0',n0,1)
3694       call readi(mcmcard,'N1',n1,6)
3695       call readi(mcmcard,'N2',n2,4)
3696       call readi(mcmcard,'N3',n3,0)
3697       call readi(mcmcard,'N4',n4,0)
3698       call readi(mcmcard,'N5',n5,0)
3699       call readi(mcmcard,'N6',n6,10)
3700       call readi(mcmcard,'N7',n7,0)
3701       call readi(mcmcard,'N8',n8,0)
3702       call readi(mcmcard,'N9',n9,0)
3703       call readi(mcmcard,'N14',n14,0)
3704       call readi(mcmcard,'N15',n15,0)
3705       call readi(mcmcard,'N16',n16,0)
3706       call readi(mcmcard,'N17',n17,0)
3707       call readi(mcmcard,'N18',n18,0)
3708
3709       vdisulf=(index(mcmcard,'DYNSS').gt.0)
3710
3711       call readi(mcmcard,'NDIFF',ndiff,2)
3712       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3713       call readi(mcmcard,'IS1',is1,1)
3714       call readi(mcmcard,'IS2',is2,8)
3715       call readi(mcmcard,'NRAN0',nran0,4)
3716       call readi(mcmcard,'NRAN1',nran1,2)
3717       call readi(mcmcard,'IRR',irr,1)
3718       call readi(mcmcard,'NSEED',nseed,20)
3719       call readi(mcmcard,'NTOTAL',ntotal,10000)
3720       call reada(mcmcard,'CUT1',cut1,2.0d0)
3721       call reada(mcmcard,'CUT2',cut2,5.0d0)
3722       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3723       call readi(mcmcard,'ICMAX',icmax,3)
3724       call readi(mcmcard,'IRESTART',irestart,0)
3725 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
3726       ntbankm=0
3727 !!bankt
3728       call reada(mcmcard,'DELE',dele,20.0d0)
3729       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3730       call readi(mcmcard,'IREF',iref,0)
3731       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3732       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3733       call readi(mcmcard,'NCONF_IN',nconf_in,0)
3734       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3735       write (iout,*) "NCONF_IN",nconf_in
3736       return
3737       end subroutine csaread
3738 !-----------------------------------------------------------------------------
3739       subroutine mcmread
3740
3741       use mcm_data
3742       use control_data, only: MaxMoveType
3743       use MD_data
3744       use minim_data
3745 !      implicit real*8 (a-h,o-z)
3746 !      include 'DIMENSIONS'
3747 !      include 'COMMON.MCM'
3748 !      include 'COMMON.MCE'
3749 !      include 'COMMON.IOUNITS'
3750 !      character(len=80) :: ucase
3751       character(len=320) :: mcmcard
3752 !el local variables
3753       integer :: i
3754 !     real(kind=8) :: var,ene
3755
3756       call card_concat(mcmcard,.true.)
3757       call readi(mcmcard,'MAXACC',maxacc,100)
3758       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3759       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3760       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3761       call readi(mcmcard,'MAXREPM',maxrepm,200)
3762       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3763       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3764       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3765       call reada(mcmcard,'E_UP',e_up,5.0D0)
3766       call reada(mcmcard,'DELTE',delte,0.1D0)
3767       call readi(mcmcard,'NSWEEP',nsweep,5)
3768       call readi(mcmcard,'NSTEPH',nsteph,0)
3769       call readi(mcmcard,'NSTEPC',nstepc,0)
3770       call reada(mcmcard,'TMIN',tmin,298.0D0)
3771       call reada(mcmcard,'TMAX',tmax,298.0D0)
3772       call readi(mcmcard,'NWINDOW',nwindow,0)
3773       call readi(mcmcard,'PRINT_MC',print_mc,0)
3774       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3775       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3776       ent_read=(index(mcmcard,'ENT_READ').gt.0)
3777       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3778       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3779       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3780       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3781       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3782       if (nwindow.gt.0) then
3783         allocate(winstart(nwindow))     !!el (maxres)
3784         allocate(winend(nwindow))       !!el
3785         allocate(winlen(nwindow))       !!el
3786         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3787         do i=1,nwindow
3788           winlen(i)=winend(i)-winstart(i)+1
3789         enddo
3790       endif
3791       if (tmax.lt.tmin) tmax=tmin
3792       if (tmax.eq.tmin) then
3793         nstepc=0
3794         nsteph=0
3795       endif
3796       if (nstepc.gt.0 .and. nsteph.gt.0) then
3797         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3798         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3799       endif
3800       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3801 ! Probabilities of different move types
3802       sumpro_type(0)=0.0D0
3803       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3804       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3805       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3806       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
3807       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3808       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3809       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3810       do i=1,MaxMoveType
3811         print *,'i',i,' sumprotype',sumpro_type(i)
3812         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3813         print *,'i',i,' sumprotype',sumpro_type(i)
3814       enddo
3815       return
3816       end subroutine mcmread
3817 !-----------------------------------------------------------------------------
3818       subroutine read_minim
3819
3820       use minim_data
3821 !      implicit real*8 (a-h,o-z)
3822 !      include 'DIMENSIONS'
3823 !      include 'COMMON.MINIM'
3824 !      include 'COMMON.IOUNITS'
3825 !      character(len=80) :: ucase
3826       character(len=320) :: minimcard
3827 !el local variables
3828 !     integer :: ntf,ik,iw_pdb
3829 !     real(kind=8) :: var,ene
3830
3831       call card_concat(minimcard,.true.)
3832       call readi(minimcard,'MAXMIN',maxmin,2000)
3833       call readi(minimcard,'MAXFUN',maxfun,5000)
3834       call readi(minimcard,'MINMIN',minmin,maxmin)
3835       call readi(minimcard,'MINFUN',minfun,maxmin)
3836       call reada(minimcard,'TOLF',tolf,1.0D-2)
3837       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3838       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3839       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3840       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3841       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3842                'Options in energy minimization:'
3843       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3844        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3845        'MinMin:',MinMin,' MinFun:',MinFun,&
3846        ' TolF:',TolF,' RTolF:',RTolF
3847       return
3848       end subroutine read_minim
3849 !-----------------------------------------------------------------------------
3850       subroutine openunits
3851
3852       use MD_data, only: usampl
3853       use csa_data
3854       use MPI_data
3855       use control_data, only:out1file
3856       use control, only: getenv_loc
3857 !      implicit real*8 (a-h,o-z)
3858 !      include 'DIMENSIONS'
3859 #ifdef MPI
3860       include 'mpif.h'
3861       character(len=16) :: form,nodename
3862       integer :: nodelen,ierror,npos
3863 #endif
3864 !      include 'COMMON.SETUP'
3865 !      include 'COMMON.IOUNITS'
3866 !      include 'COMMON.MD'
3867 !      include 'COMMON.CONTROL'
3868       integer :: lenpre,lenpot,lentmp   !,ilen
3869 !el      external ilen
3870       character(len=3) :: out1file_text !,ucase
3871       character(len=3) :: ll
3872 !el      external ucase
3873 !el local variables
3874 !     integer :: ntf,ik,iw_pdb
3875 !     real(kind=8) :: var,ene
3876 !
3877 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3878       call getenv_loc("PREFIX",prefix)
3879       pref_orig = prefix
3880       call getenv_loc("POT",pot)
3881       call getenv_loc("DIRTMP",tmpdir)
3882       call getenv_loc("CURDIR",curdir)
3883       call getenv_loc("OUT1FILE",out1file_text)
3884 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3885       out1file_text=ucase(out1file_text)
3886       if (out1file_text(1:1).eq."Y") then
3887         out1file=.true.
3888       else
3889         out1file=fg_rank.gt.0
3890       endif
3891       lenpre=ilen(prefix)
3892       lenpot=ilen(pot)
3893       lentmp=ilen(tmpdir)
3894       if (lentmp.gt.0) then
3895           write (*,'(80(1h!))')
3896           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
3897           write (*,'(80(1h!))')
3898           write (*,*)"All output files will be on node /tmp directory."
3899 #ifdef MPI
3900         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3901         if (me.eq.king) then
3902           write (*,*) "The master node is ",nodename
3903         else if (fg_rank.eq.0) then
3904           write (*,*) "I am the CG slave node ",nodename
3905         else
3906           write (*,*) "I am the FG slave node ",nodename
3907         endif
3908 #endif
3909         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3910         lenpre = lentmp+lenpre+1
3911       endif
3912       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3913 ! Get the names and open the input files
3914 #if defined(WINIFL) || defined(WINPGI)
3915       open(1,file=pref_orig(:ilen(pref_orig))// &
3916         '.inp',status='old',readonly,shared)
3917        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3918 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3919 ! Get parameter filenames and open the parameter files.
3920       call getenv_loc('BONDPAR',bondname)
3921       open (ibond,file=bondname,status='old',readonly,shared)
3922       call getenv_loc('THETPAR',thetname)
3923       open (ithep,file=thetname,status='old',readonly,shared)
3924       call getenv_loc('ROTPAR',rotname)
3925       open (irotam,file=rotname,status='old',readonly,shared)
3926       call getenv_loc('TORPAR',torname)
3927       open (itorp,file=torname,status='old',readonly,shared)
3928       call getenv_loc('TORDPAR',tordname)
3929       open (itordp,file=tordname,status='old',readonly,shared)
3930       call getenv_loc('FOURIER',fouriername)
3931       open (ifourier,file=fouriername,status='old',readonly,shared)
3932       call getenv_loc('ELEPAR',elename)
3933       open (ielep,file=elename,status='old',readonly,shared)
3934       call getenv_loc('SIDEPAR',sidename)
3935       open (isidep,file=sidename,status='old',readonly,shared)
3936 #elif (defined CRAY) || (defined AIX)
3937       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3938         action='read')
3939 !      print *,"Processor",myrank," opened file 1"
3940       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3941 !      print *,"Processor",myrank," opened file 9"
3942 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3943 ! Get parameter filenames and open the parameter files.
3944       call getenv_loc('BONDPAR',bondname)
3945       open (ibond,file=bondname,status='old',action='read')
3946 !      print *,"Processor",myrank," opened file IBOND"
3947       call getenv_loc('THETPAR',thetname)
3948       open (ithep,file=thetname,status='old',action='read')
3949 !      print *,"Processor",myrank," opened file ITHEP"
3950       call getenv_loc('ROTPAR',rotname)
3951       open (irotam,file=rotname,status='old',action='read')
3952 !      print *,"Processor",myrank," opened file IROTAM"
3953       call getenv_loc('TORPAR',torname)
3954       open (itorp,file=torname,status='old',action='read')
3955 !      print *,"Processor",myrank," opened file ITORP"
3956       call getenv_loc('TORDPAR',tordname)
3957       open (itordp,file=tordname,status='old',action='read')
3958 !      print *,"Processor",myrank," opened file ITORDP"
3959       call getenv_loc('SCCORPAR',sccorname)
3960       open (isccor,file=sccorname,status='old',action='read')
3961 !      print *,"Processor",myrank," opened file ISCCOR"
3962       call getenv_loc('FOURIER',fouriername)
3963       open (ifourier,file=fouriername,status='old',action='read')
3964 !      print *,"Processor",myrank," opened file IFOURIER"
3965       call getenv_loc('ELEPAR',elename)
3966       open (ielep,file=elename,status='old',action='read')
3967 !      print *,"Processor",myrank," opened file IELEP"
3968       call getenv_loc('SIDEPAR',sidename)
3969       open (isidep,file=sidename,status='old',action='read')
3970       call getenv_loc('LIPTRANPAR',liptranname)
3971       open (iliptranpar,file=liptranname,status='old',action='read')
3972       call getenv_loc('TUBEPAR',tubename)
3973       open (itube,file=tubename,status='old',action='read')
3974
3975 !      print *,"Processor",myrank," opened file ISIDEP"
3976 !      print *,"Processor",myrank," opened parameter files"
3977 #elif (defined G77)
3978       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3979       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3980 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3981 ! Get parameter filenames and open the parameter files.
3982       call getenv_loc('BONDPAR',bondname)
3983       open (ibond,file=bondname,status='old')
3984       call getenv_loc('THETPAR',thetname)
3985       open (ithep,file=thetname,status='old')
3986       call getenv_loc('ROTPAR',rotname)
3987       open (irotam,file=rotname,status='old')
3988       call getenv_loc('TORPAR',torname)
3989       open (itorp,file=torname,status='old')
3990       call getenv_loc('TORDPAR',tordname)
3991       open (itordp,file=tordname,status='old')
3992       call getenv_loc('SCCORPAR',sccorname)
3993       open (isccor,file=sccorname,status='old')
3994       call getenv_loc('FOURIER',fouriername)
3995       open (ifourier,file=fouriername,status='old')
3996       call getenv_loc('ELEPAR',elename)
3997       open (ielep,file=elename,status='old')
3998       call getenv_loc('SIDEPAR',sidename)
3999       open (isidep,file=sidename,status='old')
4000       call getenv_loc('LIPTRANPAR',liptranname)
4001       open (iliptranpar,file=liptranname,status='old')
4002       call getenv_loc('TUBEPAR',tubename)
4003       open (itube,file=tubename,status='old')
4004 #else
4005       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4006         readonly)
4007        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4008 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4009 ! Get parameter filenames and open the parameter files.
4010       call getenv_loc('BONDPAR',bondname)
4011       open (ibond,file=bondname,status='old',action='read')
4012       call getenv_loc('THETPAR',thetname)
4013       open (ithep,file=thetname,status='old',action='read')
4014       call getenv_loc('ROTPAR',rotname)
4015       open (irotam,file=rotname,status='old',action='read')
4016       call getenv_loc('TORPAR',torname)
4017       open (itorp,file=torname,status='old',action='read')
4018       call getenv_loc('TORDPAR',tordname)
4019       open (itordp,file=tordname,status='old',action='read')
4020       call getenv_loc('SCCORPAR',sccorname)
4021       open (isccor,file=sccorname,status='old',action='read')
4022 #ifndef CRYST_THETA
4023       call getenv_loc('THETPARPDB',thetname_pdb)
4024       print *,"thetname_pdb ",thetname_pdb
4025       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4026       print *,ithep_pdb," opened"
4027 #endif
4028       call getenv_loc('FOURIER',fouriername)
4029       open (ifourier,file=fouriername,status='old',readonly)
4030       call getenv_loc('ELEPAR',elename)
4031       open (ielep,file=elename,status='old',readonly)
4032       call getenv_loc('SIDEPAR',sidename)
4033       open (isidep,file=sidename,status='old',readonly)
4034       call getenv_loc('LIPTRANPAR',liptranname)
4035       open (iliptranpar,file=liptranname,status='old',action='read')
4036       call getenv_loc('TUBEPAR',tubename)
4037       open (itube,file=tubename,status='old',action='read')
4038
4039 #ifndef CRYST_SC
4040       call getenv_loc('ROTPARPDB',rotname_pdb)
4041       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4042 #endif
4043 #endif
4044 #ifndef OLDSCP
4045 !
4046 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4047 ! Use -DOLDSCP to use hard-coded constants instead.
4048 !
4049       call getenv_loc('SCPPAR',scpname)
4050 #if defined(WINIFL) || defined(WINPGI)
4051       open (iscpp,file=scpname,status='old',readonly,shared)
4052 #elif (defined CRAY)  || (defined AIX)
4053       open (iscpp,file=scpname,status='old',action='read')
4054 #elif (defined G77)
4055       open (iscpp,file=scpname,status='old')
4056 #else
4057       open (iscpp,file=scpname,status='old',action='read')
4058 #endif
4059 #endif
4060       call getenv_loc('PATTERN',patname)
4061 #if defined(WINIFL) || defined(WINPGI)
4062       open (icbase,file=patname,status='old',readonly,shared)
4063 #elif (defined CRAY)  || (defined AIX)
4064       open (icbase,file=patname,status='old',action='read')
4065 #elif (defined G77)
4066       open (icbase,file=patname,status='old')
4067 #else
4068       open (icbase,file=patname,status='old',action='read')
4069 #endif
4070 #ifdef MPI
4071 ! Open output file only for CG processes
4072 !      print *,"Processor",myrank," fg_rank",fg_rank
4073       if (fg_rank.eq.0) then
4074
4075       if (nodes.eq.1) then
4076         npos=3
4077       else
4078         npos = dlog10(dfloat(nodes-1))+1
4079       endif
4080       if (npos.lt.3) npos=3
4081       write (liczba,'(i1)') npos
4082       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4083         //')'
4084       write (liczba,form) me
4085       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4086         liczba(:ilen(liczba))
4087       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4088         //'.int'
4089       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4090         //'.pdb'
4091       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4092         liczba(:ilen(liczba))//'.mol2'
4093       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4094         liczba(:ilen(liczba))//'.stat'
4095       if (lentmp.gt.0) &
4096         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4097             //liczba(:ilen(liczba))//'.stat')
4098       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4099         //'.rst'
4100       if(usampl) then
4101           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4102        liczba(:ilen(liczba))//'.const'
4103       endif
4104
4105       endif
4106 #else
4107       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4108       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4109       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4110       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4111       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4112       if (lentmp.gt.0) &
4113         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4114           '.stat')
4115       rest2name=prefix(:ilen(prefix))//'.rst'
4116       if(usampl) then
4117          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4118       endif
4119 #endif
4120 #if defined(AIX) || defined(PGI)
4121       if (me.eq.king .or. .not. out1file) &
4122          open(iout,file=outname,status='unknown')
4123 #ifdef DEBUG
4124       if (fg_rank.gt.0) then
4125         write (liczba,'(i3.3)') myrank/nfgtasks
4126         write (ll,'(bz,i3.3)') fg_rank
4127         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4128          status='unknown')
4129       endif
4130 #endif
4131       if(me.eq.king) then
4132        open(igeom,file=intname,status='unknown',position='append')
4133        open(ipdb,file=pdbname,status='unknown')
4134        open(imol2,file=mol2name,status='unknown')
4135        open(istat,file=statname,status='unknown',position='append')
4136       else
4137 !1out       open(iout,file=outname,status='unknown')
4138       endif
4139 #else
4140       if (me.eq.king .or. .not.out1file) &
4141           open(iout,file=outname,status='unknown')
4142 #ifdef DEBUG
4143       if (fg_rank.gt.0) then
4144         write (liczba,'(i3.3)') myrank/nfgtasks
4145         write (ll,'(bz,i3.3)') fg_rank
4146         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4147          status='unknown')
4148       endif
4149 #endif
4150       if(me.eq.king) then
4151        open(igeom,file=intname,status='unknown',access='append')
4152        open(ipdb,file=pdbname,status='unknown')
4153        open(imol2,file=mol2name,status='unknown')
4154        open(istat,file=statname,status='unknown',access='append')
4155       else
4156 !1out       open(iout,file=outname,status='unknown')
4157       endif
4158 #endif
4159       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4160       csa_seed=prefix(:lenpre)//'.CSA.seed'
4161       csa_history=prefix(:lenpre)//'.CSA.history'
4162       csa_bank=prefix(:lenpre)//'.CSA.bank'
4163       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4164       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4165       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4166 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4167       csa_int=prefix(:lenpre)//'.int'
4168       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4169       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4170       csa_in=prefix(:lenpre)//'.CSA.in'
4171 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4172 ! Write file names
4173       if (me.eq.king)then
4174       write (iout,'(80(1h-))')
4175       write (iout,'(30x,a)') "FILE ASSIGNMENT"
4176       write (iout,'(80(1h-))')
4177       write (iout,*) "Input file                      : ",&
4178         pref_orig(:ilen(pref_orig))//'.inp'
4179       write (iout,*) "Output file                     : ",&
4180         outname(:ilen(outname))
4181       write (iout,*)
4182       write (iout,*) "Sidechain potential file        : ",&
4183         sidename(:ilen(sidename))
4184 #ifndef OLDSCP
4185       write (iout,*) "SCp potential file              : ",&
4186         scpname(:ilen(scpname))
4187 #endif
4188       write (iout,*) "Electrostatic potential file    : ",&
4189         elename(:ilen(elename))
4190       write (iout,*) "Cumulant coefficient file       : ",&
4191         fouriername(:ilen(fouriername))
4192       write (iout,*) "Torsional parameter file        : ",&
4193         torname(:ilen(torname))
4194       write (iout,*) "Double torsional parameter file : ",&
4195         tordname(:ilen(tordname))
4196       write (iout,*) "SCCOR parameter file : ",&
4197         sccorname(:ilen(sccorname))
4198       write (iout,*) "Bond & inertia constant file    : ",&
4199         bondname(:ilen(bondname))
4200       write (iout,*) "Bending parameter file          : ",&
4201         thetname(:ilen(thetname))
4202       write (iout,*) "Rotamer parameter file          : ",&
4203         rotname(:ilen(rotname))
4204 !el----
4205 #ifndef CRYST_THETA
4206       write (iout,*) "Thetpdb parameter file          : ",&
4207         thetname_pdb(:ilen(thetname_pdb))
4208 #endif
4209 !el
4210       write (iout,*) "Threading database              : ",&
4211         patname(:ilen(patname))
4212       if (lentmp.ne.0) &
4213       write (iout,*)" DIRTMP                          : ",&
4214         tmpdir(:lentmp)
4215       write (iout,'(80(1h-))')
4216       endif
4217       return
4218       end subroutine openunits
4219 !-----------------------------------------------------------------------------
4220       subroutine readrst
4221
4222       use geometry_data, only: nres,dc
4223       use MD_data
4224 !      implicit real*8 (a-h,o-z)
4225 !      include 'DIMENSIONS'
4226 !      include 'COMMON.CHAIN'
4227 !      include 'COMMON.IOUNITS'
4228 !      include 'COMMON.MD'
4229 !el local variables
4230       integer ::i,j
4231 !     real(kind=8) :: var,ene
4232
4233       open(irest2,file=rest2name,status='unknown')
4234       read(irest2,*) totT,EK,potE,totE,t_bath
4235       totTafm=totT
4236 !      do i=1,2*nres
4237 ! AL 4/17/17: Now reading d_t(0,:) too
4238       do i=0,2*nres
4239          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4240       enddo
4241 !      do i=1,2*nres
4242 ! AL 4/17/17: Now reading d_c(0,:) too
4243       do i=0,2*nres
4244          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4245       enddo
4246       if(usampl) then
4247              read (irest2,*) iset
4248       endif
4249       close(irest2)
4250       return
4251       end subroutine readrst
4252 !-----------------------------------------------------------------------------
4253       subroutine read_fragments
4254
4255       use energy_data
4256 !      use geometry
4257       use control_data, only:out1file
4258       use MD_data
4259       use MPI_data
4260 !      implicit real*8 (a-h,o-z)
4261 !      include 'DIMENSIONS'
4262 #ifdef MPI
4263       include 'mpif.h'
4264 #endif
4265 !      include 'COMMON.SETUP'
4266 !      include 'COMMON.CHAIN'
4267 !      include 'COMMON.IOUNITS'
4268 !      include 'COMMON.MD'
4269 !      include 'COMMON.CONTROL'
4270 !el local variables
4271       integer :: i
4272 !     real(kind=8) :: var,ene
4273
4274       read(inp,*) nset,nfrag,npair,nfrag_back
4275
4276 !el from module energy
4277 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
4278       if(.not.allocated(wfrag_back)) then
4279         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4280         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4281
4282         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4283         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
4284
4285         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4286         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
4287       endif
4288
4289       if(me.eq.king.or..not.out1file) &
4290        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4291         " nfrag_back",nfrag_back
4292       do iset=1,nset
4293          read(inp,*) mset(iset)
4294        do i=1,nfrag
4295          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4296            qinfrag(i,iset)
4297          if(me.eq.king.or..not.out1file) &
4298           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4299            ifrag(2,i,iset), qinfrag(i,iset)
4300        enddo
4301        do i=1,npair
4302         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4303           qinpair(i,iset)
4304         if(me.eq.king.or..not.out1file) &
4305          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4306           ipair(2,i,iset), qinpair(i,iset)
4307        enddo
4308        do i=1,nfrag_back
4309         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4310            wfrag_back(3,i,iset),&
4311            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4312         if(me.eq.king.or..not.out1file) &
4313          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4314          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4315        enddo
4316       enddo
4317       return
4318       end subroutine read_fragments
4319 !-----------------------------------------------------------------------------
4320 ! shift.F       io_csa
4321 !-----------------------------------------------------------------------------
4322       subroutine csa_read
4323
4324       use csa_data
4325 !      implicit real*8 (a-h,o-z)
4326 !      include 'DIMENSIONS'
4327 !      include 'COMMON.CSA'
4328 !      include 'COMMON.BANK'
4329 !      include 'COMMON.IOUNITS'
4330 !el local variables
4331 !     integer :: ntf,ik,iw_pdb
4332 !     real(kind=8) :: var,ene
4333
4334       open(icsa_in,file=csa_in,status="old",err=100)
4335        read(icsa_in,*) nconf
4336        read(icsa_in,*) jstart,jend
4337        read(icsa_in,*) nstmax
4338        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4339        read(icsa_in,*) nran0,nran1,irr
4340        read(icsa_in,*) nseed
4341        read(icsa_in,*) ntotal,cut1,cut2
4342        read(icsa_in,*) estop
4343        read(icsa_in,*) icmax,irestart
4344        read(icsa_in,*) ntbankm,dele,difcut
4345        read(icsa_in,*) iref,rmscut,pnccut
4346        read(icsa_in,*) ndiff
4347       close(icsa_in)
4348
4349       return
4350
4351  100  continue
4352       return
4353       end subroutine csa_read
4354 !-----------------------------------------------------------------------------
4355       subroutine initial_write
4356
4357       use csa_data
4358 !      implicit real*8 (a-h,o-z)
4359 !      include 'DIMENSIONS'
4360 !      include 'COMMON.CSA'
4361 !      include 'COMMON.BANK'
4362 !      include 'COMMON.IOUNITS'
4363 !el local variables
4364 !     integer :: ntf,ik,iw_pdb
4365 !     real(kind=8) :: var,ene
4366
4367       open(icsa_seed,file=csa_seed,status="unknown")
4368        write(icsa_seed,*) "seed"
4369       close(31)
4370 #if defined(AIX) || defined(PGI)
4371        open(icsa_history,file=csa_history,status="unknown",&
4372         position="append")
4373 #else
4374        open(icsa_history,file=csa_history,status="unknown",&
4375         access="append")
4376 #endif
4377        write(icsa_history,*) nconf
4378        write(icsa_history,*) jstart,jend
4379        write(icsa_history,*) nstmax
4380        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4381        write(icsa_history,*) nran0,nran1,irr
4382        write(icsa_history,*) nseed
4383        write(icsa_history,*) ntotal,cut1,cut2
4384        write(icsa_history,*) estop
4385        write(icsa_history,*) icmax,irestart
4386        write(icsa_history,*) ntbankm,dele,difcut
4387        write(icsa_history,*) iref,rmscut,pnccut
4388        write(icsa_history,*) ndiff
4389
4390        write(icsa_history,*)
4391       close(icsa_history)
4392
4393       open(icsa_bank1,file=csa_bank1,status="unknown")
4394        write(icsa_bank1,*) 0
4395       close(icsa_bank1)
4396
4397       return
4398       end subroutine initial_write
4399 !-----------------------------------------------------------------------------
4400       subroutine restart_write
4401
4402       use csa_data
4403 !      implicit real*8 (a-h,o-z)
4404 !      include 'DIMENSIONS'
4405 !      include 'COMMON.IOUNITS'
4406 !      include 'COMMON.CSA'
4407 !      include 'COMMON.BANK'
4408 !el local variables
4409 !     integer :: ntf,ik,iw_pdb
4410 !     real(kind=8) :: var,ene
4411
4412 #if defined(AIX) || defined(PGI)
4413        open(icsa_history,file=csa_history,position="append")
4414 #else
4415        open(icsa_history,file=csa_history,access="append")
4416 #endif
4417        write(icsa_history,*)
4418        write(icsa_history,*) "This is restart"
4419        write(icsa_history,*)
4420        write(icsa_history,*) nconf
4421        write(icsa_history,*) jstart,jend
4422        write(icsa_history,*) nstmax
4423        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4424        write(icsa_history,*) nran0,nran1,irr
4425        write(icsa_history,*) nseed
4426        write(icsa_history,*) ntotal,cut1,cut2
4427        write(icsa_history,*) estop
4428        write(icsa_history,*) icmax,irestart
4429        write(icsa_history,*) ntbankm,dele,difcut
4430        write(icsa_history,*) iref,rmscut,pnccut
4431        write(icsa_history,*) ndiff
4432        write(icsa_history,*)
4433        write(icsa_history,*) "irestart is: ", irestart
4434
4435        write(icsa_history,*)
4436       close(icsa_history)
4437
4438       return
4439       end subroutine restart_write
4440 !-----------------------------------------------------------------------------
4441 ! test.F
4442 !-----------------------------------------------------------------------------
4443       subroutine write_pdb(npdb,titelloc,ee)
4444
4445 !      implicit real*8 (a-h,o-z)
4446 !      include 'DIMENSIONS'
4447 !      include 'COMMON.IOUNITS'
4448       character(len=50) :: titelloc1
4449       character*(*) :: titelloc
4450       character(len=3) :: zahl
4451       character(len=5) :: liczba5
4452       real(kind=8) :: ee
4453       integer :: npdb   !,ilen
4454 !el      external ilen
4455 !el local variables
4456       integer :: lenpre
4457 !     real(kind=8) :: var,ene
4458
4459       titelloc1=titelloc
4460       lenpre=ilen(prefix)
4461       if (npdb.lt.1000) then
4462        call numstr(npdb,zahl)
4463        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4464       else
4465         if (npdb.lt.10000) then
4466          write(liczba5,'(i1,i4)') 0,npdb
4467         else
4468          write(liczba5,'(i5)') npdb
4469         endif
4470         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4471       endif
4472       call pdbout(ee,titelloc1,ipdb)
4473       close(ipdb)
4474       return
4475       end subroutine write_pdb
4476 !-----------------------------------------------------------------------------
4477 ! thread.F
4478 !-----------------------------------------------------------------------------
4479       subroutine write_thread_summary
4480 ! Thread the sequence through a database of known structures
4481       use control_data, only: refstr
4482 !      use geometry
4483       use energy_data, only: n_ene_comp
4484       use compare_data
4485 !      implicit real*8 (a-h,o-z)
4486 !      include 'DIMENSIONS'
4487 #ifdef MPI
4488       use MPI_data      !include 'COMMON.INFO'
4489       include 'mpif.h'
4490 #endif
4491 !      include 'COMMON.CONTROL'
4492 !      include 'COMMON.CHAIN'
4493 !      include 'COMMON.DBASE'
4494 !      include 'COMMON.INTERACT'
4495 !      include 'COMMON.VAR'
4496 !      include 'COMMON.THREAD'
4497 !      include 'COMMON.FFIELD'
4498 !      include 'COMMON.SBRIDGE'
4499 !      include 'COMMON.HEADER'
4500 !      include 'COMMON.NAMES'
4501 !      include 'COMMON.IOUNITS'
4502 !      include 'COMMON.TIME1'
4503
4504       integer,dimension(maxthread) :: ip
4505       real(kind=8),dimension(0:n_ene) :: energia
4506 !el local variables
4507       integer :: i,j,ii,jj,ipj,ik,kk,ist
4508       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4509
4510       write (iout,'(30x,a/)') &
4511        '  *********** Summary threading statistics ************'
4512       write (iout,'(a)') 'Initial energies:'
4513       write (iout,'(a4,2x,a12,14a14,3a8)') &
4514        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4515        'RMSnat','NatCONT','NNCONT','RMS'
4516 ! Energy sort patterns
4517       do i=1,nthread
4518         ip(i)=i
4519       enddo
4520       do i=1,nthread-1
4521         enet=ener(n_ene-1,ip(i))
4522         jj=i
4523         do j=i+1,nthread
4524           if (ener(n_ene-1,ip(j)).lt.enet) then
4525             jj=j
4526             enet=ener(n_ene-1,ip(j))
4527           endif
4528         enddo
4529         if (jj.ne.i) then
4530           ipj=ip(jj)
4531           ip(jj)=ip(i)
4532           ip(i)=ipj
4533         endif
4534       enddo
4535       do ik=1,nthread
4536         i=ip(ik)
4537         ii=ipatt(1,i)
4538         ist=nres_base(2,ii)+ipatt(2,i)
4539         do kk=1,n_ene_comp
4540           energia(i)=ener0(kk,i)
4541         enddo
4542         etot=ener0(n_ene_comp+1,i)
4543         rmsnat=ener0(n_ene_comp+2,i)
4544         rms=ener0(n_ene_comp+3,i)
4545         frac=ener0(n_ene_comp+4,i)
4546         frac_nn=ener0(n_ene_comp+5,i)
4547
4548         if (refstr) then
4549         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4550         i,str_nam(ii),ist+1,&
4551         (energia(print_order(kk)),kk=1,nprint_ene),&
4552         etot,rmsnat,frac,frac_nn,rms
4553         else
4554         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4555         i,str_nam(ii),ist+1,&
4556         (energia(print_order(kk)),kk=1,nprint_ene),etot
4557         endif
4558       enddo
4559       write (iout,'(//a)') 'Final energies:'
4560       write (iout,'(a4,2x,a12,17a14,3a8)') &
4561        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4562        'RMSnat','NatCONT','NNCONT','RMS'
4563       do ik=1,nthread
4564         i=ip(ik)
4565         ii=ipatt(1,i)
4566         ist=nres_base(2,ii)+ipatt(2,i)
4567         do kk=1,n_ene_comp
4568           energia(kk)=ener(kk,ik)
4569         enddo
4570         etot=ener(n_ene_comp+1,i)
4571         rmsnat=ener(n_ene_comp+2,i)
4572         rms=ener(n_ene_comp+3,i)
4573         frac=ener(n_ene_comp+4,i)
4574         frac_nn=ener(n_ene_comp+5,i)
4575         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4576         i,str_nam(ii),ist+1,&
4577         (energia(print_order(kk)),kk=1,nprint_ene),&
4578         etot,rmsnat,frac,frac_nn,rms
4579       enddo
4580       write (iout,'(/a/)') 'IEXAM array:'
4581       write (iout,'(i5)') nexcl
4582       do i=1,nexcl
4583         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4584       enddo
4585       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4586        'Max. time for threading step ',max_time_for_thread,&
4587        'Average time for threading step: ',ave_time_for_thread
4588       return
4589       end subroutine write_thread_summary
4590 !-----------------------------------------------------------------------------
4591       subroutine write_stat_thread(ithread,ipattern,ist)
4592
4593       use energy_data, only: n_ene_comp
4594       use compare_data
4595 !      implicit real*8 (a-h,o-z)
4596 !      include "DIMENSIONS"
4597 !      include "COMMON.CONTROL"
4598 !      include "COMMON.IOUNITS"
4599 !      include "COMMON.THREAD"
4600 !      include "COMMON.FFIELD"
4601 !      include "COMMON.DBASE"
4602 !      include "COMMON.NAMES"
4603       real(kind=8),dimension(0:n_ene) :: energia
4604 !el local variables
4605       integer :: ithread,ipattern,ist,i
4606       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4607
4608 #if defined(AIX) || defined(PGI)
4609       open(istat,file=statname,position='append')
4610 #else
4611       open(istat,file=statname,access='append')
4612 #endif
4613       do i=1,n_ene_comp
4614         energia(i)=ener(i,ithread)
4615       enddo
4616       etot=ener(n_ene_comp+1,ithread)
4617       rmsnat=ener(n_ene_comp+2,ithread)
4618       rms=ener(n_ene_comp+3,ithread)
4619       frac=ener(n_ene_comp+4,ithread)
4620       frac_nn=ener(n_ene_comp+5,ithread)
4621       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4622         ithread,str_nam(ipattern),ist+1,&
4623         (energia(print_order(i)),i=1,nprint_ene),&
4624         etot,rmsnat,frac,frac_nn,rms
4625       close (istat)
4626       return
4627       end subroutine write_stat_thread
4628 !-----------------------------------------------------------------------------
4629 #endif
4630 !-----------------------------------------------------------------------------
4631       end module io_config