source/unres/io_config.f90

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