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