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