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