source/unres/io_config.f90

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