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
2526 !      j=1
2527        do j=1,ntyp_molec(1) ! without U then we will take T for U
2528         write (*,*) "Im in scpho ", i, " ", j
2529         read(isidep_scpho,*) &
2530         eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2531         chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2532          write(*,*) "eps",eps_scpho(j)
2533         read(isidep_scpho,*) &
2534        (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2535        chis_scpho(j,1),chis_scpho(j,2)
2536         read(isidep_scpho,*) &
2537        (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2538         read(isidep_scpho,*) &
2539          epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j)
2540
2541        END DO
2542       allocate(aa_scpho(ntyp_molec(1)))
2543       allocate(bb_scpho(ntyp_molec(1)))
2544
2545        do j=1,ntyp_molec(1)
2546           epsij=eps_scpho(j)
2547           rrij=sigma_scpho(j)
2548 !          r0(i,j)=rrij
2549 !          r0(j,i)=rrij
2550           rrij=rrij**expon
2551 !          epsij=eps(i,j)
2552           sigeps=dsign(1.0D0,epsij)
2553           epsij=dabs(epsij)
2554           aa_scpho(j)=epsij*rrij*rrij
2555           bb_scpho(j)=-sigeps*epsij*rrij
2556         enddo
2557
2558
2559         read(isidep_peppho,*) &
2560         eps_peppho,sigma_peppho
2561         read(isidep_peppho,*) &
2562        (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2563         read(isidep_peppho,*) &
2564        (wqdip_peppho(k),k=1,2)
2565
2566           epsij=eps_peppho
2567           rrij=sigma_peppho
2568 !          r0(i,j)=rrij
2569 !          r0(j,i)=rrij
2570           rrij=rrij**expon
2571 !          epsij=eps(i,j)
2572           sigeps=dsign(1.0D0,epsij)
2573           epsij=dabs(epsij)
2574           aa_peppho=epsij*rrij*rrij
2575           bb_peppho=-sigeps*epsij*rrij
2576
2577
2578       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2579       bad(:,:)=0.0D0
2580
2581 #ifdef OLDSCP
2582 !
2583 ! Define the SC-p interaction constants (hard-coded; old style)
2584 !
2585       do i=1,ntyp
2586 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2587 ! helix formation)
2588 !       aad(i,1)=0.3D0*4.0D0**12
2589 ! Following line for constants currently implemented
2590 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2591         aad(i,1)=1.5D0*4.0D0**12
2592 !       aad(i,1)=0.17D0*5.6D0**12
2593         aad(i,2)=aad(i,1)
2594 ! "Soft" SC-p repulsion
2595         bad(i,1)=0.0D0
2596 ! Following line for constants currently implemented
2597 !       aad(i,1)=0.3D0*4.0D0**6
2598 ! "Hard" SC-p repulsion
2599         bad(i,1)=3.0D0*4.0D0**6
2600 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2601         bad(i,2)=bad(i,1)
2602 !       aad(i,1)=0.0D0
2603 !       aad(i,2)=0.0D0
2604 !       bad(i,1)=1228.8D0
2605 !       bad(i,2)=1228.8D0
2606       enddo
2607 #else
2608 !
2609 ! 8/9/01 Read the SC-p interaction constants from file
2610 !
2611       do i=1,ntyp
2612         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2613       enddo
2614       do i=1,ntyp
2615         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2616         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2617         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2618         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2619       enddo
2620 !      lprint=.true.
2621       if (lprint) then
2622         write (iout,*) "Parameters of SC-p interactions:"
2623         do i=1,ntyp
2624           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2625            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2626         enddo
2627       endif
2628 !      lprint=.false.
2629 #endif
2630       allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2631
2632       do i=1,ntyp_molec(2)
2633         read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2634       enddo
2635       do i=1,ntyp_molec(2)
2636         aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2637         bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2638       enddo
2639       r0pp=1.12246204830937298142*5.16158
2640       epspp=4.95713/4
2641       AEES=108.661
2642       BEES=0.433246
2643
2644 !
2645 ! Define the constants of the disulfide bridge
2646 !
2647       ebr=-5.50D0
2648 !
2649 ! Old arbitrary potential - commented out.
2650 !
2651 !      dbr= 4.20D0
2652 !      fbr= 3.30D0
2653 !
2654 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2655 ! energy surface of diethyl disulfide.
2656 ! A. Liwo and U. Kozlowska, 11/24/03
2657 !
2658       D0CM = 3.78d0
2659       AKCM = 15.1d0
2660       AKTH = 11.0d0
2661       AKCT = 12.0d0
2662       V1SS =-1.08d0
2663       V2SS = 7.61d0
2664       V3SS = 13.7d0
2665 !      akcm=0.0d0
2666 !      akth=0.0d0
2667 !      akct=0.0d0
2668 !      v1ss=0.0d0
2669 !      v2ss=0.0d0
2670 !      v3ss=0.0d0
2671
2672       if(me.eq.king) then
2673       write (iout,'(/a)') "Disulfide bridge parameters:"
2674       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2675       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2676       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2677       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2678         ' v3ss:',v3ss
2679       endif
2680       return
2681   111 write (iout,*) "Error reading bending energy parameters."
2682       goto 999
2683   112 write (iout,*) "Error reading rotamer energy parameters."
2684       goto 999
2685   113 write (iout,*) "Error reading torsional energy parameters."
2686       goto 999
2687   114 write (iout,*) "Error reading double torsional energy parameters."
2688       goto 999
2689   115 write (iout,*) &
2690         "Error reading cumulant (multibody energy) parameters."
2691       goto 999
2692   116 write (iout,*) "Error reading electrostatic energy parameters."
2693       goto 999
2694   117 write (iout,*) "Error reading side chain interaction parameters."
2695       goto 999
2696   118 write (iout,*) "Error reading SCp interaction parameters."
2697       goto 999
2698   119 write (iout,*) "Error reading SCCOR parameters"
2699   999 continue
2700 #ifdef MPI
2701       call MPI_Finalize(Ierror)
2702 #endif
2703       stop
2704       return
2705       end subroutine parmread
2706 #endif
2707 !-----------------------------------------------------------------------------
2708 ! printmat.f
2709 !-----------------------------------------------------------------------------
2710       subroutine printmat(ldim,m,n,iout,key,a)
2711
2712       integer :: n,ldim
2713       character(len=3),dimension(n) :: key
2714       real(kind=8),dimension(ldim,n) :: a
2715 !el local variables
2716       integer :: i,j,k,m,iout,nlim
2717
2718       do 1 i=1,n,8
2719       nlim=min0(i+7,n)
2720       write (iout,1000) (key(k),k=i,nlim)
2721       write (iout,1020)
2722  1000 format (/5x,8(6x,a3))
2723  1020 format (/80(1h-)/)
2724       do 2 j=1,n
2725       write (iout,1010) key(j),(a(j,k),k=i,nlim)
2726     2 continue
2727     1 continue
2728  1010 format (a3,2x,8(f9.4))
2729       return
2730       end subroutine printmat
2731 !-----------------------------------------------------------------------------
2732 ! readpdb.F
2733 !-----------------------------------------------------------------------------
2734       subroutine readpdb
2735 ! Read the PDB file and convert the peptide geometry into virtual-chain
2736 ! geometry.
2737       use geometry_data
2738       use energy_data, only: itype,istype
2739       use control_data
2740       use compare_data
2741       use MPI_data
2742       use control, only: rescode,sugarcode
2743 !      implicit real*8 (a-h,o-z)
2744 !      include 'DIMENSIONS'
2745 !      include 'COMMON.LOCAL'
2746 !      include 'COMMON.VAR'
2747 !      include 'COMMON.CHAIN'
2748 !      include 'COMMON.INTERACT'
2749 !      include 'COMMON.IOUNITS'
2750 !      include 'COMMON.GEO'
2751 !      include 'COMMON.NAMES'
2752 !      include 'COMMON.CONTROL'
2753 !      include 'COMMON.DISTFIT'
2754 !      include 'COMMON.SETUP'
2755       integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2756 !        ishift_pdb
2757       logical :: lprn=.true.,fail
2758       real(kind=8),dimension(3) :: e1,e2,e3
2759       real(kind=8) :: dcj,efree_temp
2760       character(len=3) :: seq,res
2761       character(len=5) :: atom
2762       character(len=80) :: card
2763       real(kind=8),dimension(3,20) :: sccor
2764       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin       !rescode,
2765       integer :: isugar,molecprev,firstion
2766       character*1 :: sugar
2767       real(kind=8) :: cou
2768       real(kind=8),dimension(3) :: ccc
2769 !el local varables
2770       integer,dimension(2,maxres/3) :: hfrag_alloc
2771       integer,dimension(4,maxres/3) :: bfrag_alloc
2772       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2773       real(kind=8),dimension(:,:), allocatable  :: c_temporary
2774       integer,dimension(:,:) , allocatable  :: itype_temporary
2775       integer,dimension(:),allocatable :: istype_temp
2776       efree_temp=0.0d0
2777       ibeg=1
2778       ishift1=0
2779       ishift=0
2780       molecule=1
2781       counter=0
2782 !      write (2,*) "UNRES_PDB",unres_pdb
2783       ires=0
2784       ires_old=0
2785       nres=0
2786       iii=0
2787       lsecondary=.false.
2788       nhfrag=0
2789       nbfrag=0
2790 !-----------------------------
2791       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2792       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2793
2794       do i=1,100000
2795         read (ipdbin,'(a80)',end=10) card
2796 !       write (iout,'(a)') card
2797         if (card(:5).eq.'HELIX') then
2798           nhfrag=nhfrag+1
2799           lsecondary=.true.
2800           read(card(22:25),*) hfrag(1,nhfrag)
2801           read(card(34:37),*) hfrag(2,nhfrag)
2802         endif
2803         if (card(:5).eq.'SHEET') then
2804           nbfrag=nbfrag+1
2805           lsecondary=.true.
2806           read(card(24:26),*) bfrag(1,nbfrag)
2807           read(card(35:37),*) bfrag(2,nbfrag)
2808 !rc----------------------------------------
2809 !rc  to be corrected !!!
2810           bfrag(3,nbfrag)=bfrag(1,nbfrag)
2811           bfrag(4,nbfrag)=bfrag(2,nbfrag)
2812 !rc----------------------------------------
2813         endif
2814         if (card(:3).eq.'END') then
2815           goto 10
2816         else if (card(:3).eq.'TER') then
2817 ! End current chain
2818           ires_old=ires+2
2819           ishift=ishift+1
2820           ishift1=ishift1+1
2821           itype(ires_old,molecule)=ntyp1_molec(molecule)
2822           itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2823           nres_molec(molecule)=nres_molec(molecule)+2
2824           ibeg=2
2825 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2826           if (unres_pdb) then
2827             do j=1,3
2828               dc(j,ires)=sccor(j,iii)
2829             enddo
2830           else
2831             call sccenter(ires,iii,sccor)
2832 !          iii=0
2833           endif
2834           iii=0
2835         endif
2836 ! Read free energy
2837         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2838 ! Fish out the ATOM cards.
2839 !        write(iout,*) 'card',card(1:20)
2840         if (index(card(1:4),'ATOM').gt.0) then
2841           read (card(12:16),*) atom
2842 !          write (iout,*) "! ",atom," !",ires
2843 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2844           read (card(23:26),*) ires
2845           read (card(18:20),'(a3)') res
2846 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2847 !     &      " ires_old",ires_old
2848 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2849 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2850           if (ires-ishift+ishift1.ne.ires_old) then
2851 ! Calculate the CM of the preceding residue.
2852 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2853             if (ibeg.eq.0) then
2854 !              write (iout,*) "Calculating sidechain center iii",iii
2855               if (unres_pdb) then
2856                 do j=1,3
2857                   dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2858                 enddo
2859               else
2860                 call sccenter(ires_old,iii,sccor)
2861               endif
2862               iii=0
2863             endif
2864 ! Start new residue.
2865             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2866               ires=ires_old
2867               cycle
2868             else if (ibeg.eq.1) then
2869               write (iout,*) "BEG ires",ires
2870               ishift=ires-1
2871               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2872                 ishift=ishift-1
2873                 itype(1,1)=ntyp1
2874                 nres_molec(molecule)=nres_molec(molecule)+1
2875               endif
2876               ires=ires-ishift+ishift1
2877               ires_old=ires
2878 !              write (iout,*) "ishift",ishift," ires",ires,&
2879 !               " ires_old",ires_old
2880               ibeg=0
2881             else if (ibeg.eq.2) then
2882 ! Start a new chain
2883               ishift=-ires_old+ires-1 !!!!!
2884               ishift1=ishift1-1    !!!!!
2885 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2886               ires=ires-ishift+ishift1
2887 !              print *,ires,ishift,ishift1
2888               ires_old=ires
2889               ibeg=0
2890             else
2891               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2892               ires=ires-ishift+ishift1
2893               ires_old=ires
2894             endif
2895 !            print *,'atom',ires,atom
2896             if (res.eq.'ACE' .or. res.eq.'NHE') then
2897               itype(ires,1)=10
2898             else
2899              if (atom.eq.'CA  '.or.atom.eq.'N   ') then
2900              molecule=1
2901               itype(ires,molecule)=rescode(ires,res,0,molecule)
2902               firstion=0
2903 !              nres_molec(molecule)=nres_molec(molecule)+1
2904             else
2905              molecule=2
2906               itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
2907 !              nres_molec(molecule)=nres_molec(molecule)+1
2908              read (card(19:19),'(a1)') sugar
2909              isugar=sugarcode(sugar,ires)
2910 !            if (ibeg.eq.1) then
2911 !              istype(1)=isugar
2912 !            else
2913               istype(ires)=isugar
2914 !              print *,"ires=",ires,istype(ires)
2915 !            endif
2916
2917             endif
2918             endif
2919           else
2920             ires=ires-ishift+ishift1
2921           endif
2922 !          write (iout,*) "ires_old",ires_old," ires",ires
2923           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2924 !            ishift1=ishift1+1
2925           endif
2926 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2927           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2928              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2929             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2930 !              print *,ires,ishift,ishift1
2931 !            write (iout,*) "backbone ",atom
2932 #ifdef DEBUG
2933             write (iout,'(2i3,2x,a,3f8.3)') &
2934             ires,itype(ires,1),res,(c(j,ires),j=1,3)
2935 #endif
2936             iii=iii+1
2937               nres_molec(molecule)=nres_molec(molecule)+1
2938             do j=1,3
2939               sccor(j,iii)=c(j,ires)
2940             enddo
2941           else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2942                atom.eq."C2'" .or. atom.eq."C3'" &
2943                .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2944             read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2945 !c            write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2946 !              print *,ires,ishift,ishift1
2947             counter=counter+1
2948 !            iii=iii+1
2949 !            do j=1,3
2950 !              sccor(j,iii)=c(j,ires)
2951 !            enddo
2952             do j=1,3
2953               c(j,ires)=c(j,ires)+ccc(j)/5.0
2954             enddo
2955              print *,counter,molecule
2956              if (counter.eq.5) then
2957 !            iii=iii+1
2958               nres_molec(molecule)=nres_molec(molecule)+1
2959               firstion=0
2960 !            do j=1,3
2961 !              sccor(j,iii)=c(j,ires)
2962 !            enddo
2963              counter=0
2964            endif
2965 !            print *, "ATOM",atom(1:3)
2966           else if (atom.eq."C5'") then
2967              read (card(19:19),'(a1)') sugar
2968              isugar=sugarcode(sugar,ires)
2969             if (ibeg.eq.1) then
2970               istype(1)=isugar
2971             else
2972               istype(ires)=isugar
2973 !              print *,ires,istype(ires)
2974             endif
2975             if (unres_pdb) then
2976               molecule=2
2977 !              print *,"nres_molec(molecule)",nres_molec(molecule),ires
2978               read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2979               nres_molec(molecule)=nres_molec(molecule)+1
2980               print *,"nres_molec(molecule)",nres_molec(molecule),ires
2981
2982             else
2983               iii=iii+1
2984               read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2985             endif
2986           else if ((atom.eq."C1'").and.unres_pdb) then
2987               iii=iii+1
2988               read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2989 !            write (*,*) card(23:27),ires,itype(ires,1)
2990           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2991                    atom.ne.'N' .and. atom.ne.'C' .and. &
2992                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2993                    atom.ne.'OXT' .and. atom(:2).ne.'3H' &
2994                    .and. atom.ne.'P  '.and. &
2995                   atom(1:1).ne.'H' .and. &
2996                   atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
2997                   ) then
2998 !            write (iout,*) "sidechain ",atom
2999 !            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3000                  if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3001 !                        write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3002             iii=iii+1
3003             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3004               endif
3005           endif
3006         else if ((ions).and.(card(1:6).eq.'HETATM')) then
3007        if (firstion.eq.0) then
3008        firstion=1
3009        if (unres_pdb) then
3010          do j=1,3
3011            dc(j,ires)=sccor(j,iii)
3012          enddo
3013        else
3014           call sccenter(ires,iii,sccor)
3015        endif
3016        endif
3017           read (card(12:16),*) atom
3018           print *,"HETATOM", atom
3019           read (card(18:20),'(a3)') res
3020           if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3021           (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
3022           .or.(atom(1:2).eq.'K ')) &
3023           then
3024            ires=ires+1
3025            if (molecule.ne.5) molecprev=molecule
3026            molecule=5
3027            nres_molec(molecule)=nres_molec(molecule)+1
3028            print *,"HERE",nres_molec(molecule)
3029            itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
3030            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3031           endif
3032         endif !atom
3033       enddo
3034    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3035       if (ires.eq.0) return
3036 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3037 ! system
3038       nres=ires
3039       if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3040         .and.(.not.unres_pdb)) &
3041          nres_molec(molecule)=nres_molec(molecule)-2
3042       print *,'I have',nres, nres_molec(:)
3043
3044       do k=1,4 ! ions are without dummy
3045        if (nres_molec(k).eq.0) cycle
3046       do i=2,nres-1
3047 !        write (iout,*) i,itype(i,1)
3048 !        if (itype(i,1).eq.ntyp1) then
3049 !          write (iout,*) "dummy",i,itype(i,1)
3050 !          do j=1,3
3051 !            c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3052 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
3053 !            dc(j,i)=c(j,i)
3054 !          enddo
3055 !        endif
3056         if (itype(i,k).eq.ntyp1_molec(k)) then
3057          if (itype(i+1,k).eq.ntyp1_molec(k)) then
3058           if (itype(i+2,k).eq.0) then
3059 !           print *,"masz sieczke"
3060            do j=1,5
3061             if (itype(i+2,j).ne.0) then
3062             itype(i+1,k)=0
3063             itype(i+1,j)=ntyp1_molec(j)
3064             nres_molec(k)=nres_molec(k)-1
3065             nres_molec(j)=nres_molec(j)+1
3066             go to 3331
3067             endif
3068            enddo
3069  3331      continue
3070           endif
3071 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3072 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3073 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3074 !           if (unres_pdb) then
3075 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3076 !            print *,i,'tu dochodze'
3077 !            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3078 !            if (fail) then
3079 !              e2(1)=0.0d0
3080 !              e2(2)=1.0d0
3081 !              e2(3)=0.0d0
3082 !            endif !fail
3083 !            print *,i,'a tu?'
3084 !            do j=1,3
3085 !             c(j,i)=c(j,i-1)-1.9d0*e2(j)
3086 !            enddo
3087 !           else   !unres_pdb
3088            do j=1,3
3089              dcj=(c(j,i-2)-c(j,i-3))/2.0
3090             if (dcj.eq.0) dcj=1.23591524223
3091              c(j,i)=c(j,i-1)+dcj
3092              c(j,nres+i)=c(j,i)
3093            enddo
3094 !          endif   !unres_pdb
3095          else     !itype(i+1,1).eq.ntyp1
3096 !          if (unres_pdb) then
3097 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3098 !            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3099 !            if (fail) then
3100 !              e2(1)=0.0d0
3101 !              e2(2)=1.0d0
3102 !              e2(3)=0.0d0
3103 !            endif
3104             do j=1,3
3105               c(j,i)=c(j,i+1)-1.9d0*e2(j)
3106             enddo
3107 !          else !unres_pdb
3108            do j=1,3
3109             dcj=(c(j,i+3)-c(j,i+2))/2.0
3110             if (dcj.eq.0) dcj=1.23591524223
3111             c(j,i)=c(j,i+1)-dcj
3112             c(j,nres+i)=c(j,i)
3113            enddo
3114 !          endif !unres_pdb
3115          endif !itype(i+1,1).eq.ntyp1
3116         endif  !itype.eq.ntyp1
3117
3118       enddo
3119       enddo
3120 ! Calculate the CM of the last side chain.
3121       if (iii.gt.0)  then
3122       if (unres_pdb) then
3123         do j=1,3
3124           dc(j,ires)=sccor(j,iii)
3125         enddo
3126       else
3127         call sccenter(ires,iii,sccor)
3128       endif
3129       endif
3130 !      nres=ires
3131       nsup=nres
3132       nstart_sup=1
3133 !      print *,"molecule",molecule
3134       if ((itype(nres,1).ne.10)) then
3135         nres=nres+1
3136           if (molecule.eq.5) molecule=molecprev
3137         itype(nres,molecule)=ntyp1_molec(molecule)
3138         nres_molec(molecule)=nres_molec(molecule)+1
3139 !        if (unres_pdb) then
3140 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3141 !          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3142 !          if (fail) then
3143 !            e2(1)=0.0d0
3144 !            e2(2)=1.0d0
3145 !            e2(3)=0.0d0
3146 !          endif
3147 !          do j=1,3
3148 !            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3149 !          enddo
3150 !        else
3151         do j=1,3
3152           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3153           c(j,nres)=c(j,nres-1)+dcj
3154           c(j,2*nres)=c(j,nres)
3155         enddo
3156 !        endif
3157       endif
3158 !     print *,'I have',nres, nres_molec(:)
3159
3160 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3161 #ifdef WHAM_RUN
3162       if (nres.ne.nres0) then
3163         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3164                        " NRES0=",nres0
3165         stop "Error nres value in WHAM input"
3166       endif
3167 #endif
3168 !---------------------------------
3169 !el reallocate tables
3170 !      do i=1,maxres/3
3171 !       do j=1,2
3172 !         hfrag_alloc(j,i)=hfrag(j,i)
3173 !        enddo
3174 !       do j=1,4
3175 !         bfrag_alloc(j,i)=bfrag(j,i)
3176 !        enddo
3177 !      enddo
3178
3179 !      deallocate(hfrag)
3180 !      deallocate(bfrag)
3181 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
3182 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3183 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3184 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
3185
3186 !      do i=1,nhfrag
3187 !       do j=1,2
3188 !         hfrag(j,i)=hfrag_alloc(j,i)
3189 !        enddo
3190 !      enddo
3191 !      do i=1,nbfrag
3192 !       do j=1,4
3193 !         bfrag(j,i)=bfrag_alloc(j,i)
3194 !        enddo
3195 !      enddo
3196 !el end reallocate tables
3197 !---------------------------------
3198       do i=2,nres-1
3199         do j=1,3
3200           c(j,i+nres)=dc(j,i)
3201         enddo
3202       enddo
3203       do j=1,3
3204         c(j,nres+1)=c(j,1)
3205         c(j,2*nres)=c(j,nres)
3206       enddo
3207
3208       if (itype(1,1).eq.ntyp1) then
3209         nsup=nsup-1
3210         nstart_sup=2
3211         if (unres_pdb) then
3212 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3213           call refsys(2,3,4,e1,e2,e3,fail)
3214           if (fail) then
3215             e2(1)=0.0d0
3216             e2(2)=1.0d0
3217             e2(3)=0.0d0
3218           endif
3219           do j=1,3
3220             c(j,1)=c(j,2)-1.9d0*e2(j)
3221           enddo
3222         else
3223         do j=1,3
3224           dcj=(c(j,4)-c(j,3))/2.0
3225           c(j,1)=c(j,2)-dcj
3226           c(j,nres+1)=c(j,1)
3227         enddo
3228         endif
3229       endif
3230 ! First lets assign correct dummy to correct type of chain
3231 ! 1) First residue
3232       if (itype(1,1).eq.ntyp1) then
3233         if (itype(2,1).eq.0) then
3234          do j=2,5
3235            if (itype(2,j).ne.0) then
3236            itype(1,1)=0
3237            itype(1,j)=ntyp1_molec(j)
3238            nres_molec(1)=nres_molec(1)-1
3239            nres_molec(j)=nres_molec(j)+1
3240            go to 3231
3241            endif
3242          enddo
3243 3231    continue
3244         endif
3245        endif
3246        print *,'I have',nres, nres_molec(:)
3247
3248 ! Copy the coordinates to reference coordinates
3249 !      do i=1,2*nres
3250 !        do j=1,3
3251 !          cref(j,i)=c(j,i)
3252 !        enddo
3253 !      enddo
3254 ! Calculate internal coordinates.
3255       if (lprn) then
3256       write (iout,'(/a)') &
3257         "Cartesian coordinates of the reference structure"
3258       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3259        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3260       do ires=1,nres
3261         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3262           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3263           (c(j,ires+nres),j=1,3)
3264       enddo
3265       endif
3266 ! znamy już nres wiec mozna alokowac tablice
3267 ! Calculate internal coordinates.
3268       if(me.eq.king.or..not.out1file)then
3269        write (iout,'(a)') &
3270          "Backbone and SC coordinates as read from the PDB"
3271        do ires=1,nres
3272         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3273           ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3274           (c(j,nres+ires),j=1,3)
3275        enddo
3276       endif
3277 ! NOW LETS ROCK! SORTING
3278       allocate(c_temporary(3,2*nres))
3279       allocate(itype_temporary(nres,5))
3280       allocate(molnum(nres+1))
3281       allocate(istype_temp(nres))
3282        itype_temporary(:,:)=0
3283       seqalingbegin=1
3284       do k=1,5
3285         do i=1,nres
3286          if (itype(i,k).ne.0) then
3287           do j=1,3
3288           c_temporary(j,seqalingbegin)=c(j,i)
3289           c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3290
3291           enddo
3292           itype_temporary(seqalingbegin,k)=itype(i,k)
3293           print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3294           istype_temp(seqalingbegin)=istype(i)
3295           molnum(seqalingbegin)=k
3296           seqalingbegin=seqalingbegin+1
3297          endif
3298         enddo
3299        enddo
3300        do i=1,2*nres
3301         do j=1,3
3302         c(j,i)=c_temporary(j,i)
3303         enddo
3304        enddo
3305        do k=1,5
3306         do i=1,nres
3307          itype(i,k)=itype_temporary(i,k)
3308          istype(i)=istype_temp(i)
3309         enddo
3310        enddo
3311       if (itype(1,1).eq.ntyp1) then
3312         nsup=nsup-1
3313         nstart_sup=2
3314         if (unres_pdb) then
3315 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3316           call refsys(2,3,4,e1,e2,e3,fail)
3317           if (fail) then
3318             e2(1)=0.0d0
3319             e2(2)=1.0d0
3320             e2(3)=0.0d0
3321           endif
3322           do j=1,3
3323             c(j,1)=c(j,2)-1.9d0*e2(j)
3324           enddo
3325         else
3326         do j=1,3
3327           dcj=(c(j,4)-c(j,3))/2.0
3328           c(j,1)=c(j,2)-dcj
3329           c(j,nres+1)=c(j,1)
3330         enddo
3331         endif
3332       endif
3333
3334       if (lprn) then
3335       write (iout,'(/a)') &
3336         "Cartesian coordinates of the reference structure after sorting"
3337       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3338        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3339       do ires=1,nres
3340         write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3341           (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3342           (c(j,ires+nres),j=1,3)
3343       enddo
3344       endif
3345
3346 !       print *,seqalingbegin,nres
3347       if(.not.allocated(vbld)) then
3348        allocate(vbld(2*nres))
3349        do i=1,2*nres
3350          vbld(i)=0.d0
3351        enddo
3352       endif
3353       if(.not.allocated(vbld_inv)) then
3354        allocate(vbld_inv(2*nres))
3355        do i=1,2*nres
3356          vbld_inv(i)=0.d0
3357        enddo
3358       endif
3359 !!!el
3360       if(.not.allocated(theta)) then
3361         allocate(theta(nres+2))
3362         theta(:)=0.0d0
3363       endif
3364
3365       if(.not.allocated(phi)) allocate(phi(nres+2))
3366       if(.not.allocated(alph)) allocate(alph(nres+2))
3367       if(.not.allocated(omeg)) allocate(omeg(nres+2))
3368       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3369       if(.not.allocated(phiref)) allocate(phiref(nres+2))
3370       if(.not.allocated(costtab)) allocate(costtab(nres))
3371       if(.not.allocated(sinttab)) allocate(sinttab(nres))
3372       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3373       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3374       if(.not.allocated(xxref)) allocate(xxref(nres))
3375       if(.not.allocated(yyref)) allocate(yyref(nres))
3376       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3377       if(.not.allocated(dc_norm)) then
3378 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3379         allocate(dc_norm(3,0:2*nres+2))
3380         dc_norm(:,:)=0.d0
3381       endif
3382
3383       call int_from_cart(.true.,.false.)
3384       call sc_loc_geom(.false.)
3385       do i=1,nres
3386         thetaref(i)=theta(i)
3387         phiref(i)=phi(i)
3388       enddo
3389 !      do i=1,2*nres
3390 !        vbld_inv(i)=0.d0
3391 !        vbld(i)=0.d0
3392 !      enddo
3393
3394       do i=1,nres-1
3395         do j=1,3
3396           dc(j,i)=c(j,i+1)-c(j,i)
3397           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3398         enddo
3399       enddo
3400       do i=2,nres-1
3401         do j=1,3
3402           dc(j,i+nres)=c(j,i+nres)-c(j,i)
3403           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3404         enddo
3405 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3406 !        vbld_inv(i+nres)
3407       enddo
3408 !      call chainbuild
3409 ! Copy the coordinates to reference coordinates
3410 ! Splits to single chain if occurs
3411
3412 !      do i=1,2*nres
3413 !        do j=1,3
3414 !          cref(j,i,cou)=c(j,i)
3415 !        enddo
3416 !      enddo
3417 !
3418       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3419       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3420       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3421 !-----------------------------
3422       kkk=1
3423       lll=0
3424       cou=1
3425         write (iout,*) "symetr", symetr
3426       do i=1,nres
3427       lll=lll+1
3428 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3429       if (i.gt.1) then
3430       if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3431       chain_length=lll-1
3432       kkk=kkk+1
3433 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3434       lll=1
3435       endif
3436       endif
3437         do j=1,3
3438           cref(j,i,cou)=c(j,i)
3439           cref(j,i+nres,cou)=c(j,i+nres)
3440           if (i.le.nres) then
3441           chain_rep(j,lll,kkk)=c(j,i)
3442           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3443           endif
3444          enddo
3445       enddo
3446       write (iout,*) chain_length
3447       if (chain_length.eq.0) chain_length=nres
3448       do j=1,3
3449       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3450       chain_rep(j,chain_length+nres,symetr) &
3451       =chain_rep(j,chain_length+nres,1)
3452       enddo
3453 ! diagnostic
3454 !       write (iout,*) "spraw lancuchy",chain_length,symetr
3455 !       do i=1,4
3456 !         do kkk=1,chain_length
3457 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3458 !         enddo
3459 !        enddo
3460 ! enddiagnostic
3461 ! makes copy of chains
3462         write (iout,*) "symetr", symetr
3463       do j=1,3
3464       dc(j,0)=c(j,1)
3465       enddo
3466
3467       if (symetr.gt.1) then
3468        call permut(symetr)
3469        nperm=1
3470        do i=1,symetr
3471        nperm=nperm*i
3472        enddo
3473        do i=1,nperm
3474        write(iout,*) (tabperm(i,kkk),kkk=1,4)
3475        enddo
3476        do i=1,nperm
3477         cou=0
3478         do kkk=1,symetr
3479          icha=tabperm(i,kkk)
3480 !         write (iout,*) i,icha
3481          do lll=1,chain_length
3482           cou=cou+1
3483            if (cou.le.nres) then
3484            do j=1,3
3485             kupa=mod(lll,chain_length)
3486             iprzes=(kkk-1)*chain_length+lll
3487             if (kupa.eq.0) kupa=chain_length
3488 !            write (iout,*) "kupa", kupa
3489             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3490             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3491           enddo
3492           endif
3493          enddo
3494         enddo
3495        enddo
3496        endif
3497 !-koniec robienia kopii
3498 ! diag
3499       do kkk=1,nperm
3500       write (iout,*) "nowa struktura", nperm
3501       do i=1,nres
3502         write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3503       cref(2,i,kkk),&
3504       cref(3,i,kkk),cref(1,nres+i,kkk),&
3505       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3506       enddo
3507   100 format (//'              alpha-carbon coordinates       ',&
3508                 '     centroid coordinates'/ &
3509                 '       ', 6X,'X',11X,'Y',11X,'Z', &
3510                                 10X,'X',11X,'Y',11X,'Z')
3511   110 format (a,'(',i3,')',6f12.5)
3512
3513       enddo
3514 !c enddiag
3515       do j=1,nbfrag
3516         do i=1,4
3517          bfrag(i,j)=bfrag(i,j)-ishift
3518         enddo
3519       enddo
3520
3521       do j=1,nhfrag
3522         do i=1,2
3523          hfrag(i,j)=hfrag(i,j)-ishift
3524         enddo
3525       enddo
3526       ishift_pdb=ishift
3527
3528       return
3529       end subroutine readpdb
3530 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3531 !-----------------------------------------------------------------------------
3532 ! readrtns_CSA.F
3533 !-----------------------------------------------------------------------------
3534       subroutine read_control
3535 !
3536 ! Read contorl data
3537 !
3538 !      use geometry_data
3539       use comm_machsw
3540       use energy_data
3541       use control_data
3542       use compare_data
3543       use MCM_data
3544       use map_data
3545       use csa_data
3546       use MD_data
3547       use MPI_data
3548       use random, only: random_init
3549 !      implicit real*8 (a-h,o-z)
3550 !      include 'DIMENSIONS'
3551 #ifdef MP
3552       use prng, only:prng_restart
3553       include 'mpif.h'
3554       logical :: OKRandom!, prng_restart
3555       real(kind=8) :: r1
3556 #endif
3557 !      include 'COMMON.IOUNITS'
3558 !      include 'COMMON.TIME1'
3559 !      include 'COMMON.THREAD'
3560 !      include 'COMMON.SBRIDGE'
3561 !      include 'COMMON.CONTROL'
3562 !      include 'COMMON.MCM'
3563 !      include 'COMMON.MAP'
3564 !      include 'COMMON.HEADER'
3565 !      include 'COMMON.CSA'
3566 !      include 'COMMON.CHAIN'
3567 !      include 'COMMON.MUCA'
3568 !      include 'COMMON.MD'
3569 !      include 'COMMON.FFIELD'
3570 !      include 'COMMON.INTERACT'
3571 !      include 'COMMON.SETUP'
3572 !el      integer :: KDIAG,ICORFL,IXDR
3573 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3574       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3575         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
3576 !      character(len=80) :: ucase
3577       character(len=640) :: controlcard
3578
3579       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3580       integer i
3581
3582       nglob_csa=0
3583       eglob_csa=1d99
3584       nmin_csa=0
3585       read (INP,'(a)') titel
3586       call card_concat(controlcard,.true.)
3587 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3588 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3589       call reada(controlcard,'SEED',seed,0.0D0)
3590       call random_init(seed)
3591 ! Set up the time limit (caution! The time must be input in minutes!)
3592       read_cart=index(controlcard,'READ_CART').gt.0
3593       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3594       call readi(controlcard,'SYM',symetr,1)
3595       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3596       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3597       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3598       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3599       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3600       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3601       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3602       call reada(controlcard,'DRMS',drms,0.1D0)
3603       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3604        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3605        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3606        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3607        write (iout,'(a,f10.1)')'DRMS    = ',drms
3608        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3609        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3610       endif
3611       call readi(controlcard,'NZ_START',nz_start,0)
3612       call readi(controlcard,'NZ_END',nz_end,0)
3613       call readi(controlcard,'IZ_SC',iz_sc,0)
3614       timlim=60.0D0*timlim
3615       safety = 60.0d0*safety
3616       timem=timlim
3617       modecalc=0
3618       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3619 !C SHIELD keyword sets if the shielding effect of side-chains is used
3620 !C 0 denots no shielding is used all peptide are equally despite the
3621 !C solvent accesible area
3622 !C 1 the newly introduced function
3623 !C 2 reseved for further possible developement
3624       call readi(controlcard,'SHIELD',shield_mode,0)
3625 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3626         write(iout,*) "shield_mode",shield_mode
3627 !C  Varibles set size of box
3628       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3629       protein=index(controlcard,"PROTEIN").gt.0
3630       ions=index(controlcard,"IONS").gt.0
3631       nucleic=index(controlcard,"NUCLEIC").gt.0
3632       write (iout,*) "with_theta_constr ",with_theta_constr
3633       AFMlog=(index(controlcard,'AFM'))
3634       selfguide=(index(controlcard,'SELFGUIDE'))
3635       print *,'AFMlog',AFMlog,selfguide,"KUPA"
3636       call readi(controlcard,'GENCONSTR',genconstr,0)
3637       call reada(controlcard,'BOXX',boxxsize,100.0d0)
3638       call reada(controlcard,'BOXY',boxysize,100.0d0)
3639       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3640       call readi(controlcard,'TUBEMOD',tubemode,0)
3641       write (iout,*) TUBEmode,"TUBEMODE"
3642       if (TUBEmode.gt.0) then
3643        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3644        call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3645        call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3646        call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3647        call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3648        call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3649        call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3650        buftubebot=bordtubebot+tubebufthick
3651        buftubetop=bordtubetop-tubebufthick
3652       endif
3653
3654 ! CUTOFFF ON ELECTROSTATICS
3655       call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3656       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3657       write(iout,*) "R_CUT_ELE=",r_cut_ele
3658 ! Lipidic parameters
3659       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3660       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3661       if (lipthick.gt.0.0d0) then
3662        bordliptop=(boxzsize+lipthick)/2.0
3663        bordlipbot=bordliptop-lipthick
3664       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3665       write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3666       buflipbot=bordlipbot+lipbufthick
3667       bufliptop=bordliptop-lipbufthick
3668       if ((lipbufthick*2.0d0).gt.lipthick) &
3669        write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3670       endif !lipthick.gt.0
3671       write(iout,*) "bordliptop=",bordliptop
3672       write(iout,*) "bordlipbot=",bordlipbot
3673       write(iout,*) "bufliptop=",bufliptop
3674       write(iout,*) "buflipbot=",buflipbot
3675       write (iout,*) "SHIELD MODE",shield_mode
3676
3677 !C-------------------------
3678       minim=(index(controlcard,'MINIMIZE').gt.0)
3679       dccart=(index(controlcard,'CART').gt.0)
3680       overlapsc=(index(controlcard,'OVERLAP').gt.0)
3681       overlapsc=.not.overlapsc
3682       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3683       searchsc=.not.searchsc
3684       sideadd=(index(controlcard,'SIDEADD').gt.0)
3685       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3686       outpdb=(index(controlcard,'PDBOUT').gt.0)
3687       outmol2=(index(controlcard,'MOL2OUT').gt.0)
3688       pdbref=(index(controlcard,'PDBREF').gt.0)
3689       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3690       indpdb=index(controlcard,'PDBSTART')
3691       extconf=(index(controlcard,'EXTCONF').gt.0)
3692       call readi(controlcard,'IPRINT',iprint,0)
3693       call readi(controlcard,'MAXGEN',maxgen,10000)
3694       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3695       call readi(controlcard,"KDIAG",kdiag,0)
3696       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3697       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3698        write (iout,*) "RESCALE_MODE",rescale_mode
3699       split_ene=index(controlcard,'SPLIT_ENE').gt.0
3700       if (index(controlcard,'REGULAR').gt.0.0D0) then
3701         call reada(controlcard,'WEIDIS',weidis,0.1D0)
3702         modecalc=1
3703         refstr=.true.
3704       endif
3705       if (index(controlcard,'CHECKGRAD').gt.0) then
3706         modecalc=5
3707         if (index(controlcard,'CART').gt.0) then
3708           icheckgrad=1
3709         elseif (index(controlcard,'CARINT').gt.0) then
3710           icheckgrad=2
3711         else
3712           icheckgrad=3
3713         endif
3714       elseif (index(controlcard,'THREAD').gt.0) then
3715         modecalc=2
3716         call readi(controlcard,'THREAD',nthread,0)
3717         if (nthread.gt.0) then
3718           call reada(controlcard,'WEIDIS',weidis,0.1D0)
3719         else
3720           if (fg_rank.eq.0) &
3721           write (iout,'(a)')'A number has to follow the THREAD keyword.'
3722           stop 'Error termination in Read_Control.'
3723         endif
3724       else if (index(controlcard,'MCMA').gt.0) then
3725         modecalc=3
3726       else if (index(controlcard,'MCEE').gt.0) then
3727         modecalc=6
3728       else if (index(controlcard,'MULTCONF').gt.0) then
3729         modecalc=4
3730       else if (index(controlcard,'MAP').gt.0) then
3731         modecalc=7
3732         call readi(controlcard,'MAP',nmap,0)
3733       else if (index(controlcard,'CSA').gt.0) then
3734         modecalc=8
3735 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3736 !rc
3737 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3738 !rc
3739 !rc        modecalc=9
3740 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3741 !fmc        modecalc=10
3742       else if (index(controlcard,'SOFTREG').gt.0) then
3743         modecalc=11
3744       else if (index(controlcard,'CHECK_BOND').gt.0) then
3745         modecalc=-1
3746       else if (index(controlcard,'TEST').gt.0) then
3747         modecalc=-2
3748       else if (index(controlcard,'MD').gt.0) then
3749         modecalc=12
3750       else if (index(controlcard,'RE ').gt.0) then
3751         modecalc=14
3752       endif
3753
3754       lmuca=index(controlcard,'MUCA').gt.0
3755       call readi(controlcard,'MUCADYN',mucadyn,0)
3756       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3757       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3758        then
3759        write (iout,*) 'MUCADYN=',mucadyn
3760        write (iout,*) 'MUCASMOOTH=',muca_smooth
3761       endif
3762
3763       iscode=index(controlcard,'ONE_LETTER')
3764       indphi=index(controlcard,'PHI')
3765       indback=index(controlcard,'BACK')
3766       iranconf=index(controlcard,'RAND_CONF')
3767       i2ndstr=index(controlcard,'USE_SEC_PRED')
3768       gradout=index(controlcard,'GRADOUT').gt.0
3769       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3770       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3771       if (me.eq.king .or. .not.out1file ) &
3772         write (iout,*) "DISTCHAINMAX",distchainmax
3773
3774       if(me.eq.king.or..not.out1file) &
3775        write (iout,'(2a)') diagmeth(kdiag),&
3776         ' routine used to diagonalize matrices.'
3777       if (shield_mode.gt.0) then
3778       pi=3.141592d0
3779 !C VSolvSphere the volume of solving sphere
3780 !C      print *,pi,"pi"
3781 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3782 !C there will be no distinction between proline peptide group and normal peptide
3783 !C group in case of shielding parameters
3784       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3785       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3786       write (iout,*) VSolvSphere,VSolvSphere_div
3787 !C long axis of side chain
3788       do i=1,ntyp
3789       long_r_sidechain(i)=vbldsc0(1,i)
3790       short_r_sidechain(i)=sigma0(i)
3791       write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3792          sigma0(i)
3793       enddo
3794       buff_shield=1.0d0
3795       endif
3796       return
3797       end subroutine read_control
3798 !-----------------------------------------------------------------------------
3799       subroutine read_REMDpar
3800 !
3801 ! Read REMD settings
3802 !
3803 !       use control
3804 !       use energy
3805 !       use geometry
3806       use REMD_data
3807       use MPI_data
3808       use control_data, only:out1file
3809 !      implicit real*8 (a-h,o-z)
3810 !      include 'DIMENSIONS'
3811 !      include 'COMMON.IOUNITS'
3812 !      include 'COMMON.TIME1'
3813 !      include 'COMMON.MD'
3814       use MD_data
3815 !el #ifndef LANG0
3816 !el      include 'COMMON.LANGEVIN'
3817 !el #else
3818 !el      include 'COMMON.LANGEVIN.lang0'
3819 !el #endif
3820 !      include 'COMMON.INTERACT'
3821 !      include 'COMMON.NAMES'
3822 !      include 'COMMON.GEO'
3823 !      include 'COMMON.REMD'
3824 !      include 'COMMON.CONTROL'
3825 !      include 'COMMON.SETUP'
3826 !      character(len=80) :: ucase
3827       character(len=320) :: controlcard
3828       character(len=3200) :: controlcard1
3829       integer :: iremd_m_total
3830 !el local variables
3831       integer :: i
3832 !     real(kind=8) :: var,ene
3833
3834       if(me.eq.king.or..not.out1file) &
3835        write (iout,*) "REMD setup"
3836
3837       call card_concat(controlcard,.true.)
3838       call readi(controlcard,"NREP",nrep,3)
3839       call readi(controlcard,"NSTEX",nstex,1000)
3840       call reada(controlcard,"RETMIN",retmin,10.0d0)
3841       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3842       mremdsync=(index(controlcard,'SYNC').gt.0)
3843       call readi(controlcard,"NSYN",i_sync_step,100)
3844       restart1file=(index(controlcard,'REST1FILE').gt.0)
3845       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3846       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3847       if(max_cache_traj_use.gt.max_cache_traj) &
3848                 max_cache_traj_use=max_cache_traj
3849       if(me.eq.king.or..not.out1file) then
3850 !d       if (traj1file) then
3851 !rc caching is in testing - NTWX is not ignored
3852 !d        write (iout,*) "NTWX value is ignored"
3853 !d        write (iout,*) "  trajectory is stored to one file by master"
3854 !d        write (iout,*) "  before exchange at NSTEX intervals"
3855 !d       endif
3856        write (iout,*) "NREP= ",nrep
3857        write (iout,*) "NSTEX= ",nstex
3858        write (iout,*) "SYNC= ",mremdsync
3859        write (iout,*) "NSYN= ",i_sync_step
3860        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3861       endif
3862       remd_tlist=.false.
3863       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3864       if (index(controlcard,'TLIST').gt.0) then
3865          remd_tlist=.true.
3866          call card_concat(controlcard1,.true.)
3867          read(controlcard1,*) (remd_t(i),i=1,nrep)
3868          if(me.eq.king.or..not.out1file) &
3869           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3870       endif
3871       remd_mlist=.false.
3872       if (index(controlcard,'MLIST').gt.0) then
3873          remd_mlist=.true.
3874          call card_concat(controlcard1,.true.)
3875          read(controlcard1,*) (remd_m(i),i=1,nrep)
3876          if(me.eq.king.or..not.out1file) then
3877           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3878           iremd_m_total=0
3879           do i=1,nrep
3880            iremd_m_total=iremd_m_total+remd_m(i)
3881           enddo
3882           write (iout,*) 'Total number of replicas ',iremd_m_total
3883          endif
3884       endif
3885       if(me.eq.king.or..not.out1file) &
3886        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3887       return
3888       end subroutine read_REMDpar
3889 !-----------------------------------------------------------------------------
3890       subroutine read_MDpar
3891 !
3892 ! Read MD settings
3893 !
3894       use control_data, only: r_cut,rlamb,out1file
3895       use energy_data
3896       use geometry_data, only: pi
3897       use MPI_data
3898 !      implicit real*8 (a-h,o-z)
3899 !      include 'DIMENSIONS'
3900 !      include 'COMMON.IOUNITS'
3901 !      include 'COMMON.TIME1'
3902 !      include 'COMMON.MD'
3903       use MD_data
3904 !el #ifndef LANG0
3905 !el      include 'COMMON.LANGEVIN'
3906 !el #else
3907 !el      include 'COMMON.LANGEVIN.lang0'
3908 !el #endif
3909 !      include 'COMMON.INTERACT'
3910 !      include 'COMMON.NAMES'
3911 !      include 'COMMON.GEO'
3912 !      include 'COMMON.SETUP'
3913 !      include 'COMMON.CONTROL'
3914 !      include 'COMMON.SPLITELE'
3915 !      character(len=80) :: ucase
3916       character(len=320) :: controlcard
3917 !el local variables
3918       integer :: i,j
3919       real(kind=8) :: eta
3920
3921       call card_concat(controlcard,.true.)
3922       call readi(controlcard,"NSTEP",n_timestep,1000000)
3923       call readi(controlcard,"NTWE",ntwe,100)
3924       call readi(controlcard,"NTWX",ntwx,1000)
3925       call reada(controlcard,"DT",d_time,1.0d-1)
3926       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3927       call reada(controlcard,"DAMAX",damax,1.0d1)
3928       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3929       call readi(controlcard,"LANG",lang,0)
3930       RESPA = index(controlcard,"RESPA") .gt. 0
3931       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3932       ntime_split0=ntime_split
3933       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3934       ntime_split0=ntime_split
3935       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3936       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3937       rest = index(controlcard,"REST").gt.0
3938       tbf = index(controlcard,"TBF").gt.0
3939       usampl = index(controlcard,"USAMPL").gt.0
3940       mdpdb = index(controlcard,"MDPDB").gt.0
3941       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3942       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3943       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3944       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3945       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3946       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3947       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3948       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3949       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3950       large = index(controlcard,"LARGE").gt.0
3951       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3952       rattle = index(controlcard,"RATTLE").gt.0
3953 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3954       nset=0
3955       if(usampl) then
3956         call read_fragments
3957       endif
3958
3959       if(me.eq.king.or..not.out1file) then
3960        write (iout,*)
3961        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3962        write (iout,*)
3963        write (iout,'(a)') "The units are:"
3964        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3965        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3966         " acceleration: angstrom/(48.9 fs)**2"
3967        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3968        write (iout,*)
3969        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3970        write (iout,'(a60,f10.5,a)') &
3971         "Initial time step of numerical integration:",d_time,&
3972         " natural units"
3973        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3974        if (RESPA) then
3975         write (iout,'(2a,i4,a)') &
3976           "A-MTS algorithm used; initial time step for fast-varying",&
3977           " short-range forces split into",ntime_split," steps."
3978         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3979          r_cut," lambda",rlamb
3980        endif
3981        write (iout,'(2a,f10.5)') &
3982         "Maximum acceleration threshold to reduce the time step",&
3983         "/increase split number:",damax
3984        write (iout,'(2a,f10.5)') &
3985         "Maximum predicted energy drift to reduce the timestep",&
3986         "/increase split number:",edriftmax
3987        write (iout,'(a60,f10.5)') &
3988        "Maximum velocity threshold to reduce velocities:",dvmax
3989        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3990        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3991        if (rattle) write (iout,'(a60)') &
3992         "Rattle algorithm used to constrain the virtual bonds"
3993       endif
3994       reset_fricmat=1000
3995       if (lang.gt.0) then
3996         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3997         call reada(controlcard,"RWAT",rwat,1.4d0)
3998         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3999         surfarea=index(controlcard,"SURFAREA").gt.0
4000         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4001         if(me.eq.king.or..not.out1file)then
4002          write (iout,'(/a,$)') "Langevin dynamics calculation"
4003          if (lang.eq.1) then
4004           write (iout,'(a/)') &
4005             " with direct integration of Langevin equations"
4006          else if (lang.eq.2) then
4007           write (iout,'(a/)') " with TINKER stochasic MD integrator"
4008          else if (lang.eq.3) then
4009           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4010          else if (lang.eq.4) then
4011           write (iout,'(a/)') " in overdamped mode"
4012          else
4013           write (iout,'(//a,i5)') &
4014             "=========== ERROR: Unknown Langevin dynamics mode:",lang
4015           stop
4016          endif
4017          write (iout,'(a60,f10.5)') "Temperature:",t_bath
4018          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4019          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4020          write (iout,'(a60,f10.5)') &
4021          "Scaling factor of the friction forces:",scal_fric
4022          if (surfarea) write (iout,'(2a,i10,a)') &
4023            "Friction coefficients will be scaled by solvent-accessible",&
4024            " surface area every",reset_fricmat," steps."
4025         endif
4026 ! Calculate friction coefficients and bounds of stochastic forces
4027         eta=6*pi*cPoise*etawat
4028         if(me.eq.king.or..not.out1file) &
4029          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4030           eta
4031 !        allocate(gamp
4032         do j=1,5 !types of molecules
4033         gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4034         stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4035         enddo
4036         allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4037         do j=1,5 !types of molecules
4038         do i=1,ntyp
4039           gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4040           stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4041         enddo
4042         enddo
4043
4044         if(me.eq.king.or..not.out1file)then
4045          write (iout,'(/2a/)') &
4046          "Radii of site types and friction coefficients and std's of",&
4047          " stochastic forces of fully exposed sites"
4048          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4049          do i=1,ntyp
4050           write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4051            gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4052          enddo
4053         endif
4054       else if (tbf) then
4055         if(me.eq.king.or..not.out1file)then
4056          write (iout,'(a)') "Berendsen bath calculation"
4057          write (iout,'(a60,f10.5)') "Temperature:",t_bath
4058          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4059          if (reset_moment) &
4060          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4061          count_reset_moment," steps"
4062          if (reset_vel) &
4063           write (iout,'(a,i10,a)') &
4064           "Velocities will be reset at random every",count_reset_vel,&
4065          " steps"
4066         endif
4067       else
4068         if(me.eq.king.or..not.out1file) &
4069          write (iout,'(a31)') "Microcanonical mode calculation"
4070       endif
4071       if(me.eq.king.or..not.out1file)then
4072        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4073        if (usampl) then
4074           write(iout,*) "MD running with constraints."
4075           write(iout,*) "Equilibration time ", eq_time, " mtus."
4076           write(iout,*) "Constraining ", nfrag," fragments."
4077           write(iout,*) "Length of each fragment, weight and q0:"
4078           do iset=1,nset
4079            write (iout,*) "Set of restraints #",iset
4080            do i=1,nfrag
4081               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4082                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4083            enddo
4084            write(iout,*) "constraints between ", npair, "fragments."
4085            write(iout,*) "constraint pairs, weights and q0:"
4086            do i=1,npair
4087             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4088                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4089            enddo
4090            write(iout,*) "angle constraints within ", nfrag_back,&
4091             "backbone fragments."
4092            write(iout,*) "fragment, weights:"
4093            do i=1,nfrag_back
4094             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4095                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4096                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4097            enddo
4098           enddo
4099         iset=mod(kolor,nset)+1
4100        endif
4101       endif
4102       if(me.eq.king.or..not.out1file) &
4103        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4104       return
4105       end subroutine read_MDpar
4106 !-----------------------------------------------------------------------------
4107       subroutine map_read
4108
4109       use map_data
4110 !      implicit real*8 (a-h,o-z)
4111 !      include 'DIMENSIONS'
4112 !      include 'COMMON.MAP'
4113 !      include 'COMMON.IOUNITS'
4114       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4115       character(len=80) :: mapcard      !,ucase
4116 !el local variables
4117       integer :: imap
4118 !     real(kind=8) :: var,ene
4119
4120       do imap=1,nmap
4121         read (inp,'(a)') mapcard
4122         mapcard=ucase(mapcard)
4123         if (index(mapcard,'PHI').gt.0) then
4124           kang(imap)=1
4125         else if (index(mapcard,'THE').gt.0) then
4126           kang(imap)=2
4127         else if (index(mapcard,'ALP').gt.0) then
4128           kang(imap)=3
4129         else if (index(mapcard,'OME').gt.0) then
4130           kang(imap)=4
4131         else
4132           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4133           stop 'Error - illegal variable spec in MAP card.'
4134         endif
4135         call readi (mapcard,'RES1',res1(imap),0)
4136         call readi (mapcard,'RES2',res2(imap),0)
4137         if (res1(imap).eq.0) then
4138           res1(imap)=res2(imap)
4139         else if (res2(imap).eq.0) then
4140           res2(imap)=res1(imap)
4141         endif
4142         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4143           write (iout,'(a)') &
4144           'Error - illegal definition of variable group in MAP.'
4145           stop 'Error - illegal definition of variable group in MAP.'
4146         endif
4147         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4148         call reada(mapcard,'TO',ang_to(imap),0.0D0)
4149         call readi(mapcard,'NSTEP',nstep(imap),0)
4150         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4151           write (iout,'(a)') &
4152            'Illegal boundary and/or step size specification in MAP.'
4153           stop 'Illegal boundary and/or step size specification in MAP.'
4154         endif
4155       enddo ! imap
4156       return
4157       end subroutine map_read
4158 !-----------------------------------------------------------------------------
4159       subroutine csaread
4160
4161       use control_data, only: vdisulf
4162       use csa_data
4163 !      implicit real*8 (a-h,o-z)
4164 !      include 'DIMENSIONS'
4165 !      include 'COMMON.IOUNITS'
4166 !      include 'COMMON.GEO'
4167 !      include 'COMMON.CSA'
4168 !      include 'COMMON.BANK'
4169 !      include 'COMMON.CONTROL'
4170 !      character(len=80) :: ucase
4171       character(len=620) :: mcmcard
4172 !el local variables
4173 !     integer :: ntf,ik,iw_pdb
4174 !     real(kind=8) :: var,ene
4175
4176       call card_concat(mcmcard,.true.)
4177
4178       call readi(mcmcard,'NCONF',nconf,50)
4179       call readi(mcmcard,'NADD',nadd,0)
4180       call readi(mcmcard,'JSTART',jstart,1)
4181       call readi(mcmcard,'JEND',jend,1)
4182       call readi(mcmcard,'NSTMAX',nstmax,500000)
4183       call readi(mcmcard,'N0',n0,1)
4184       call readi(mcmcard,'N1',n1,6)
4185       call readi(mcmcard,'N2',n2,4)
4186       call readi(mcmcard,'N3',n3,0)
4187       call readi(mcmcard,'N4',n4,0)
4188       call readi(mcmcard,'N5',n5,0)
4189       call readi(mcmcard,'N6',n6,10)
4190       call readi(mcmcard,'N7',n7,0)
4191       call readi(mcmcard,'N8',n8,0)
4192       call readi(mcmcard,'N9',n9,0)
4193       call readi(mcmcard,'N14',n14,0)
4194       call readi(mcmcard,'N15',n15,0)
4195       call readi(mcmcard,'N16',n16,0)
4196       call readi(mcmcard,'N17',n17,0)
4197       call readi(mcmcard,'N18',n18,0)
4198
4199       vdisulf=(index(mcmcard,'DYNSS').gt.0)
4200
4201       call readi(mcmcard,'NDIFF',ndiff,2)
4202       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4203       call readi(mcmcard,'IS1',is1,1)
4204       call readi(mcmcard,'IS2',is2,8)
4205       call readi(mcmcard,'NRAN0',nran0,4)
4206       call readi(mcmcard,'NRAN1',nran1,2)
4207       call readi(mcmcard,'IRR',irr,1)
4208       call readi(mcmcard,'NSEED',nseed,20)
4209       call readi(mcmcard,'NTOTAL',ntotal,10000)
4210       call reada(mcmcard,'CUT1',cut1,2.0d0)
4211       call reada(mcmcard,'CUT2',cut2,5.0d0)
4212       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4213       call readi(mcmcard,'ICMAX',icmax,3)
4214       call readi(mcmcard,'IRESTART',irestart,0)
4215 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
4216       ntbankm=0
4217 !!bankt
4218       call reada(mcmcard,'DELE',dele,20.0d0)
4219       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4220       call readi(mcmcard,'IREF',iref,0)
4221       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4222       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4223       call readi(mcmcard,'NCONF_IN',nconf_in,0)
4224       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4225       write (iout,*) "NCONF_IN",nconf_in
4226       return
4227       end subroutine csaread
4228 !-----------------------------------------------------------------------------
4229       subroutine mcmread
4230
4231       use mcm_data
4232       use control_data, only: MaxMoveType
4233       use MD_data
4234       use minim_data
4235 !      implicit real*8 (a-h,o-z)
4236 !      include 'DIMENSIONS'
4237 !      include 'COMMON.MCM'
4238 !      include 'COMMON.MCE'
4239 !      include 'COMMON.IOUNITS'
4240 !      character(len=80) :: ucase
4241       character(len=320) :: mcmcard
4242 !el local variables
4243       integer :: i
4244 !     real(kind=8) :: var,ene
4245
4246       call card_concat(mcmcard,.true.)
4247       call readi(mcmcard,'MAXACC',maxacc,100)
4248       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4249       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4250       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4251       call readi(mcmcard,'MAXREPM',maxrepm,200)
4252       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4253       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4254       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4255       call reada(mcmcard,'E_UP',e_up,5.0D0)
4256       call reada(mcmcard,'DELTE',delte,0.1D0)
4257       call readi(mcmcard,'NSWEEP',nsweep,5)
4258       call readi(mcmcard,'NSTEPH',nsteph,0)
4259       call readi(mcmcard,'NSTEPC',nstepc,0)
4260       call reada(mcmcard,'TMIN',tmin,298.0D0)
4261       call reada(mcmcard,'TMAX',tmax,298.0D0)
4262       call readi(mcmcard,'NWINDOW',nwindow,0)
4263       call readi(mcmcard,'PRINT_MC',print_mc,0)
4264       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4265       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4266       ent_read=(index(mcmcard,'ENT_READ').gt.0)
4267       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4268       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4269       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4270       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4271       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4272       if (nwindow.gt.0) then
4273         allocate(winstart(nwindow))     !!el (maxres)
4274         allocate(winend(nwindow))       !!el
4275         allocate(winlen(nwindow))       !!el
4276         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4277         do i=1,nwindow
4278           winlen(i)=winend(i)-winstart(i)+1
4279         enddo
4280       endif
4281       if (tmax.lt.tmin) tmax=tmin
4282       if (tmax.eq.tmin) then
4283         nstepc=0
4284         nsteph=0
4285       endif
4286       if (nstepc.gt.0 .and. nsteph.gt.0) then
4287         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4288         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4289       endif
4290       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4291 ! Probabilities of different move types
4292       sumpro_type(0)=0.0D0
4293       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4294       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4295       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4296       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
4297       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4298       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4299       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4300       do i=1,MaxMoveType
4301         print *,'i',i,' sumprotype',sumpro_type(i)
4302         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4303         print *,'i',i,' sumprotype',sumpro_type(i)
4304       enddo
4305       return
4306       end subroutine mcmread
4307 !-----------------------------------------------------------------------------
4308       subroutine read_minim
4309
4310       use minim_data
4311 !      implicit real*8 (a-h,o-z)
4312 !      include 'DIMENSIONS'
4313 !      include 'COMMON.MINIM'
4314 !      include 'COMMON.IOUNITS'
4315 !      character(len=80) :: ucase
4316       character(len=320) :: minimcard
4317 !el local variables
4318 !     integer :: ntf,ik,iw_pdb
4319 !     real(kind=8) :: var,ene
4320
4321       call card_concat(minimcard,.true.)
4322       call readi(minimcard,'MAXMIN',maxmin,2000)
4323       call readi(minimcard,'MAXFUN',maxfun,5000)
4324       call readi(minimcard,'MINMIN',minmin,maxmin)
4325       call readi(minimcard,'MINFUN',minfun,maxmin)
4326       call reada(minimcard,'TOLF',tolf,1.0D-2)
4327       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4328       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4329       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4330       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4331       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4332                'Options in energy minimization:'
4333       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4334        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4335        'MinMin:',MinMin,' MinFun:',MinFun,&
4336        ' TolF:',TolF,' RTolF:',RTolF
4337       return
4338       end subroutine read_minim
4339 !-----------------------------------------------------------------------------
4340       subroutine openunits
4341
4342       use MD_data, only: usampl
4343       use csa_data
4344       use MPI_data
4345       use control_data, only:out1file
4346       use control, only: getenv_loc
4347 !      implicit real*8 (a-h,o-z)
4348 !      include 'DIMENSIONS'
4349 #ifdef MPI
4350       include 'mpif.h'
4351       character(len=16) :: form,nodename
4352       integer :: nodelen,ierror,npos
4353 #endif
4354 !      include 'COMMON.SETUP'
4355 !      include 'COMMON.IOUNITS'
4356 !      include 'COMMON.MD'
4357 !      include 'COMMON.CONTROL'
4358       integer :: lenpre,lenpot,lentmp   !,ilen
4359 !el      external ilen
4360       character(len=3) :: out1file_text !,ucase
4361       character(len=3) :: ll
4362 !el      external ucase
4363 !el local variables
4364 !     integer :: ntf,ik,iw_pdb
4365 !     real(kind=8) :: var,ene
4366 !
4367 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4368       call getenv_loc("PREFIX",prefix)
4369       pref_orig = prefix
4370       call getenv_loc("POT",pot)
4371       call getenv_loc("DIRTMP",tmpdir)
4372       call getenv_loc("CURDIR",curdir)
4373       call getenv_loc("OUT1FILE",out1file_text)
4374 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4375       out1file_text=ucase(out1file_text)
4376       if (out1file_text(1:1).eq."Y") then
4377         out1file=.true.
4378       else
4379         out1file=fg_rank.gt.0
4380       endif
4381       lenpre=ilen(prefix)
4382       lenpot=ilen(pot)
4383       lentmp=ilen(tmpdir)
4384       if (lentmp.gt.0) then
4385           write (*,'(80(1h!))')
4386           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
4387           write (*,'(80(1h!))')
4388           write (*,*)"All output files will be on node /tmp directory."
4389 #ifdef MPI
4390         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4391         if (me.eq.king) then
4392           write (*,*) "The master node is ",nodename
4393         else if (fg_rank.eq.0) then
4394           write (*,*) "I am the CG slave node ",nodename
4395         else
4396           write (*,*) "I am the FG slave node ",nodename
4397         endif
4398 #endif
4399         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4400         lenpre = lentmp+lenpre+1
4401       endif
4402       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4403 ! Get the names and open the input files
4404 #if defined(WINIFL) || defined(WINPGI)
4405       open(1,file=pref_orig(:ilen(pref_orig))// &
4406         '.inp',status='old',readonly,shared)
4407        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4408 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4409 ! Get parameter filenames and open the parameter files.
4410       call getenv_loc('BONDPAR',bondname)
4411       open (ibond,file=bondname,status='old',readonly,shared)
4412       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4413       open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4414       call getenv_loc('THETPAR',thetname)
4415       open (ithep,file=thetname,status='old',readonly,shared)
4416       call getenv_loc('ROTPAR',rotname)
4417       open (irotam,file=rotname,status='old',readonly,shared)
4418       call getenv_loc('TORPAR',torname)
4419       open (itorp,file=torname,status='old',readonly,shared)
4420       call getenv_loc('TORDPAR',tordname)
4421       open (itordp,file=tordname,status='old',readonly,shared)
4422       call getenv_loc('FOURIER',fouriername)
4423       open (ifourier,file=fouriername,status='old',readonly,shared)
4424       call getenv_loc('ELEPAR',elename)
4425       open (ielep,file=elename,status='old',readonly,shared)
4426       call getenv_loc('SIDEPAR',sidename)
4427       open (isidep,file=sidename,status='old',readonly,shared)
4428
4429       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4430       open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4431       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4432       open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4433       call getenv_loc('TORPAR_NUCL',torname_nucl)
4434       open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4435       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4436       open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4437       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4438       open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4439
4440
4441 #elif (defined CRAY) || (defined AIX)
4442       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4443         action='read')
4444 !      print *,"Processor",myrank," opened file 1"
4445       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4446 !      print *,"Processor",myrank," opened file 9"
4447 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4448 ! Get parameter filenames and open the parameter files.
4449       call getenv_loc('BONDPAR',bondname)
4450       open (ibond,file=bondname,status='old',action='read')
4451       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4452       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4453
4454 !      print *,"Processor",myrank," opened file IBOND"
4455       call getenv_loc('THETPAR',thetname)
4456       open (ithep,file=thetname,status='old',action='read')
4457 !      print *,"Processor",myrank," opened file ITHEP"
4458       call getenv_loc('ROTPAR',rotname)
4459       open (irotam,file=rotname,status='old',action='read')
4460 !      print *,"Processor",myrank," opened file IROTAM"
4461       call getenv_loc('TORPAR',torname)
4462       open (itorp,file=torname,status='old',action='read')
4463 !      print *,"Processor",myrank," opened file ITORP"
4464       call getenv_loc('TORDPAR',tordname)
4465       open (itordp,file=tordname,status='old',action='read')
4466 !      print *,"Processor",myrank," opened file ITORDP"
4467       call getenv_loc('SCCORPAR',sccorname)
4468       open (isccor,file=sccorname,status='old',action='read')
4469 !      print *,"Processor",myrank," opened file ISCCOR"
4470       call getenv_loc('FOURIER',fouriername)
4471       open (ifourier,file=fouriername,status='old',action='read')
4472 !      print *,"Processor",myrank," opened file IFOURIER"
4473       call getenv_loc('ELEPAR',elename)
4474       open (ielep,file=elename,status='old',action='read')
4475 !      print *,"Processor",myrank," opened file IELEP"
4476       call getenv_loc('SIDEPAR',sidename)
4477       open (isidep,file=sidename,status='old',action='read')
4478
4479       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4480       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4481       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4482       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4483       call getenv_loc('TORPAR_NUCL',torname_nucl)
4484       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4485       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4486       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4487       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4488       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4489
4490       call getenv_loc('LIPTRANPAR',liptranname)
4491       open (iliptranpar,file=liptranname,status='old',action='read')
4492       call getenv_loc('TUBEPAR',tubename)
4493       open (itube,file=tubename,status='old',action='read')
4494       call getenv_loc('IONPAR',ionname)
4495       open (iion,file=ionname,status='old',action='read')
4496
4497 !      print *,"Processor",myrank," opened file ISIDEP"
4498 !      print *,"Processor",myrank," opened parameter files"
4499 #elif (defined G77)
4500       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4501       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4502 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4503 ! Get parameter filenames and open the parameter files.
4504       call getenv_loc('BONDPAR',bondname)
4505       open (ibond,file=bondname,status='old')
4506       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4507       open (ibond_nucl,file=bondname_nucl,status='old')
4508
4509       call getenv_loc('THETPAR',thetname)
4510       open (ithep,file=thetname,status='old')
4511       call getenv_loc('ROTPAR',rotname)
4512       open (irotam,file=rotname,status='old')
4513       call getenv_loc('TORPAR',torname)
4514       open (itorp,file=torname,status='old')
4515       call getenv_loc('TORDPAR',tordname)
4516       open (itordp,file=tordname,status='old')
4517       call getenv_loc('SCCORPAR',sccorname)
4518       open (isccor,file=sccorname,status='old')
4519       call getenv_loc('FOURIER',fouriername)
4520       open (ifourier,file=fouriername,status='old')
4521       call getenv_loc('ELEPAR',elename)
4522       open (ielep,file=elename,status='old')
4523       call getenv_loc('SIDEPAR',sidename)
4524       open (isidep,file=sidename,status='old')
4525
4526       open (ithep_nucl,file=thetname_nucl,status='old')
4527       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4528       open (irotam_nucl,file=rotname_nucl,status='old')
4529       call getenv_loc('TORPAR_NUCL',torname_nucl)
4530       open (itorp_nucl,file=torname_nucl,status='old')
4531       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4532       open (itordp_nucl,file=tordname_nucl,status='old')
4533       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4534       open (isidep_nucl,file=sidename_nucl,status='old')
4535
4536       call getenv_loc('LIPTRANPAR',liptranname)
4537       open (iliptranpar,file=liptranname,status='old')
4538       call getenv_loc('TUBEPAR',tubename)
4539       open (itube,file=tubename,status='old')
4540       call getenv_loc('IONPAR',ionname)
4541       open (iion,file=ionname,status='old')
4542 #else
4543       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4544         readonly)
4545        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4546 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4547 ! Get parameter filenames and open the parameter files.
4548       call getenv_loc('BONDPAR',bondname)
4549       open (ibond,file=bondname,status='old',action='read')
4550       call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4551       open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4552       call getenv_loc('THETPAR',thetname)
4553       open (ithep,file=thetname,status='old',action='read')
4554       call getenv_loc('ROTPAR',rotname)
4555       open (irotam,file=rotname,status='old',action='read')
4556       call getenv_loc('TORPAR',torname)
4557       open (itorp,file=torname,status='old',action='read')
4558       call getenv_loc('TORDPAR',tordname)
4559       open (itordp,file=tordname,status='old',action='read')
4560       call getenv_loc('SCCORPAR',sccorname)
4561       open (isccor,file=sccorname,status='old',action='read')
4562 #ifndef CRYST_THETA
4563       call getenv_loc('THETPARPDB',thetname_pdb)
4564       print *,"thetname_pdb ",thetname_pdb
4565       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4566       print *,ithep_pdb," opened"
4567 #endif
4568       call getenv_loc('FOURIER',fouriername)
4569       open (ifourier,file=fouriername,status='old',readonly)
4570       call getenv_loc('ELEPAR',elename)
4571       open (ielep,file=elename,status='old',readonly)
4572       call getenv_loc('SIDEPAR',sidename)
4573       open (isidep,file=sidename,status='old',readonly)
4574
4575       call getenv_loc('THETPAR_NUCL',thetname_nucl)
4576       open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4577       call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4578       open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4579       call getenv_loc('TORPAR_NUCL',torname_nucl)
4580       open (itorp_nucl,file=torname_nucl,status='old',action='read')
4581       call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4582       open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4583       call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4584       open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4585       call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4586       open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4587       call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4588       open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4589       call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4590       open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4591       call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4592       open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4593
4594
4595       call getenv_loc('LIPTRANPAR',liptranname)
4596       open (iliptranpar,file=liptranname,status='old',action='read')
4597       call getenv_loc('TUBEPAR',tubename)
4598       open (itube,file=tubename,status='old',action='read')
4599       call getenv_loc('IONPAR',ionname)
4600       open (iion,file=ionname,status='old',action='read')
4601
4602 #ifndef CRYST_SC
4603       call getenv_loc('ROTPARPDB',rotname_pdb)
4604       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4605 #endif
4606 #endif
4607       call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4608 #if defined(WINIFL) || defined(WINPGI)
4609       open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4610 #elif (defined CRAY)  || (defined AIX)
4611       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4612 #elif (defined G77)
4613       open (iscpp_nucl,file=scpname_nucl,status='old')
4614 #else
4615       open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4616 #endif
4617
4618 #ifndef OLDSCP
4619 !
4620 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4621 ! Use -DOLDSCP to use hard-coded constants instead.
4622 !
4623       call getenv_loc('SCPPAR',scpname)
4624 #if defined(WINIFL) || defined(WINPGI)
4625       open (iscpp,file=scpname,status='old',readonly,shared)
4626 #elif (defined CRAY)  || (defined AIX)
4627       open (iscpp,file=scpname,status='old',action='read')
4628 #elif (defined G77)
4629       open (iscpp,file=scpname,status='old')
4630 #else
4631       open (iscpp,file=scpname,status='old',action='read')
4632 #endif
4633 #endif
4634       call getenv_loc('PATTERN',patname)
4635 #if defined(WINIFL) || defined(WINPGI)
4636       open (icbase,file=patname,status='old',readonly,shared)
4637 #elif (defined CRAY)  || (defined AIX)
4638       open (icbase,file=patname,status='old',action='read')
4639 #elif (defined G77)
4640       open (icbase,file=patname,status='old')
4641 #else
4642       open (icbase,file=patname,status='old',action='read')
4643 #endif
4644 #ifdef MPI
4645 ! Open output file only for CG processes
4646 !      print *,"Processor",myrank," fg_rank",fg_rank
4647       if (fg_rank.eq.0) then
4648
4649       if (nodes.eq.1) then
4650         npos=3
4651       else
4652         npos = dlog10(dfloat(nodes-1))+1
4653       endif
4654       if (npos.lt.3) npos=3
4655       write (liczba,'(i1)') npos
4656       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4657         //')'
4658       write (liczba,form) me
4659       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4660         liczba(:ilen(liczba))
4661       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4662         //'.int'
4663       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4664         //'.pdb'
4665       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4666         liczba(:ilen(liczba))//'.mol2'
4667       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4668         liczba(:ilen(liczba))//'.stat'
4669       if (lentmp.gt.0) &
4670         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4671             //liczba(:ilen(liczba))//'.stat')
4672       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4673         //'.rst'
4674       if(usampl) then
4675           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4676        liczba(:ilen(liczba))//'.const'
4677       endif
4678
4679       endif
4680 #else
4681       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4682       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4683       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4684       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4685       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4686       if (lentmp.gt.0) &
4687         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4688           '.stat')
4689       rest2name=prefix(:ilen(prefix))//'.rst'
4690       if(usampl) then
4691          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4692       endif
4693 #endif
4694 #if defined(AIX) || defined(PGI)
4695       if (me.eq.king .or. .not. out1file) &
4696          open(iout,file=outname,status='unknown')
4697 #ifdef DEBUG
4698       if (fg_rank.gt.0) then
4699         write (liczba,'(i3.3)') myrank/nfgtasks
4700         write (ll,'(bz,i3.3)') fg_rank
4701         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4702          status='unknown')
4703       endif
4704 #endif
4705       if(me.eq.king) then
4706        open(igeom,file=intname,status='unknown',position='append')
4707        open(ipdb,file=pdbname,status='unknown')
4708        open(imol2,file=mol2name,status='unknown')
4709        open(istat,file=statname,status='unknown',position='append')
4710       else
4711 !1out       open(iout,file=outname,status='unknown')
4712       endif
4713 #else
4714       if (me.eq.king .or. .not.out1file) &
4715           open(iout,file=outname,status='unknown')
4716 #ifdef DEBUG
4717       if (fg_rank.gt.0) then
4718         write (liczba,'(i3.3)') myrank/nfgtasks
4719         write (ll,'(bz,i3.3)') fg_rank
4720         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4721          status='unknown')
4722       endif
4723 #endif
4724       if(me.eq.king) then
4725        open(igeom,file=intname,status='unknown',access='append')
4726        open(ipdb,file=pdbname,status='unknown')
4727        open(imol2,file=mol2name,status='unknown')
4728        open(istat,file=statname,status='unknown',access='append')
4729       else
4730 !1out       open(iout,file=outname,status='unknown')
4731       endif
4732 #endif
4733       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4734       csa_seed=prefix(:lenpre)//'.CSA.seed'
4735       csa_history=prefix(:lenpre)//'.CSA.history'
4736       csa_bank=prefix(:lenpre)//'.CSA.bank'
4737       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4738       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4739       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4740 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4741       csa_int=prefix(:lenpre)//'.int'
4742       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4743       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4744       csa_in=prefix(:lenpre)//'.CSA.in'
4745 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4746 ! Write file names
4747       if (me.eq.king)then
4748       write (iout,'(80(1h-))')
4749       write (iout,'(30x,a)') "FILE ASSIGNMENT"
4750       write (iout,'(80(1h-))')
4751       write (iout,*) "Input file                      : ",&
4752         pref_orig(:ilen(pref_orig))//'.inp'
4753       write (iout,*) "Output file                     : ",&
4754         outname(:ilen(outname))
4755       write (iout,*)
4756       write (iout,*) "Sidechain potential file        : ",&
4757         sidename(:ilen(sidename))
4758 #ifndef OLDSCP
4759       write (iout,*) "SCp potential file              : ",&
4760         scpname(:ilen(scpname))
4761 #endif
4762       write (iout,*) "Electrostatic potential file    : ",&
4763         elename(:ilen(elename))
4764       write (iout,*) "Cumulant coefficient file       : ",&
4765         fouriername(:ilen(fouriername))
4766       write (iout,*) "Torsional parameter file        : ",&
4767         torname(:ilen(torname))
4768       write (iout,*) "Double torsional parameter file : ",&
4769         tordname(:ilen(tordname))
4770       write (iout,*) "SCCOR parameter file : ",&
4771         sccorname(:ilen(sccorname))
4772       write (iout,*) "Bond & inertia constant file    : ",&
4773         bondname(:ilen(bondname))
4774       write (iout,*) "Bending parameter file          : ",&
4775         thetname(:ilen(thetname))
4776       write (iout,*) "Rotamer parameter file          : ",&
4777         rotname(:ilen(rotname))
4778 !el----
4779 #ifndef CRYST_THETA
4780       write (iout,*) "Thetpdb parameter file          : ",&
4781         thetname_pdb(:ilen(thetname_pdb))
4782 #endif
4783 !el
4784       write (iout,*) "Threading database              : ",&
4785         patname(:ilen(patname))
4786       if (lentmp.ne.0) &
4787       write (iout,*)" DIRTMP                          : ",&
4788         tmpdir(:lentmp)
4789       write (iout,'(80(1h-))')
4790       endif
4791       return
4792       end subroutine openunits
4793 !-----------------------------------------------------------------------------
4794       subroutine readrst
4795
4796       use geometry_data, only: nres,dc
4797       use MD_data
4798 !      implicit real*8 (a-h,o-z)
4799 !      include 'DIMENSIONS'
4800 !      include 'COMMON.CHAIN'
4801 !      include 'COMMON.IOUNITS'
4802 !      include 'COMMON.MD'
4803 !el local variables
4804       integer ::i,j
4805 !     real(kind=8) :: var,ene
4806
4807       open(irest2,file=rest2name,status='unknown')
4808       read(irest2,*) totT,EK,potE,totE,t_bath
4809       totTafm=totT
4810 !      do i=1,2*nres
4811 ! AL 4/17/17: Now reading d_t(0,:) too
4812       do i=0,2*nres
4813          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4814       enddo
4815 !      do i=1,2*nres
4816 ! AL 4/17/17: Now reading d_c(0,:) too
4817       do i=0,2*nres
4818          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4819       enddo
4820       if(usampl) then
4821              read (irest2,*) iset
4822       endif
4823       close(irest2)
4824       return
4825       end subroutine readrst
4826 !-----------------------------------------------------------------------------
4827       subroutine read_fragments
4828
4829       use energy_data
4830 !      use geometry
4831       use control_data, only:out1file
4832       use MD_data
4833       use MPI_data
4834 !      implicit real*8 (a-h,o-z)
4835 !      include 'DIMENSIONS'
4836 #ifdef MPI
4837       include 'mpif.h'
4838 #endif
4839 !      include 'COMMON.SETUP'
4840 !      include 'COMMON.CHAIN'
4841 !      include 'COMMON.IOUNITS'
4842 !      include 'COMMON.MD'
4843 !      include 'COMMON.CONTROL'
4844 !el local variables
4845       integer :: i
4846 !     real(kind=8) :: var,ene
4847
4848       read(inp,*) nset,nfrag,npair,nfrag_back
4849
4850 !el from module energy
4851 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
4852       if(.not.allocated(wfrag_back)) then
4853         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4854         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4855
4856         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4857         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
4858
4859         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4860         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
4861       endif
4862
4863       if(me.eq.king.or..not.out1file) &
4864        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4865         " nfrag_back",nfrag_back
4866       do iset=1,nset
4867          read(inp,*) mset(iset)
4868        do i=1,nfrag
4869          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4870            qinfrag(i,iset)
4871          if(me.eq.king.or..not.out1file) &
4872           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4873            ifrag(2,i,iset), qinfrag(i,iset)
4874        enddo
4875        do i=1,npair
4876         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4877           qinpair(i,iset)
4878         if(me.eq.king.or..not.out1file) &
4879          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4880           ipair(2,i,iset), qinpair(i,iset)
4881        enddo
4882        do i=1,nfrag_back
4883         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4884            wfrag_back(3,i,iset),&
4885            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4886         if(me.eq.king.or..not.out1file) &
4887          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4888          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4889        enddo
4890       enddo
4891       return
4892       end subroutine read_fragments
4893 !-----------------------------------------------------------------------------
4894 ! shift.F       io_csa
4895 !-----------------------------------------------------------------------------
4896       subroutine csa_read
4897
4898       use csa_data
4899 !      implicit real*8 (a-h,o-z)
4900 !      include 'DIMENSIONS'
4901 !      include 'COMMON.CSA'
4902 !      include 'COMMON.BANK'
4903 !      include 'COMMON.IOUNITS'
4904 !el local variables
4905 !     integer :: ntf,ik,iw_pdb
4906 !     real(kind=8) :: var,ene
4907
4908       open(icsa_in,file=csa_in,status="old",err=100)
4909        read(icsa_in,*) nconf
4910        read(icsa_in,*) jstart,jend
4911        read(icsa_in,*) nstmax
4912        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4913        read(icsa_in,*) nran0,nran1,irr
4914        read(icsa_in,*) nseed
4915        read(icsa_in,*) ntotal,cut1,cut2
4916        read(icsa_in,*) estop
4917        read(icsa_in,*) icmax,irestart
4918        read(icsa_in,*) ntbankm,dele,difcut
4919        read(icsa_in,*) iref,rmscut,pnccut
4920        read(icsa_in,*) ndiff
4921       close(icsa_in)
4922
4923       return
4924
4925  100  continue
4926       return
4927       end subroutine csa_read
4928 !-----------------------------------------------------------------------------
4929       subroutine initial_write
4930
4931       use csa_data
4932 !      implicit real*8 (a-h,o-z)
4933 !      include 'DIMENSIONS'
4934 !      include 'COMMON.CSA'
4935 !      include 'COMMON.BANK'
4936 !      include 'COMMON.IOUNITS'
4937 !el local variables
4938 !     integer :: ntf,ik,iw_pdb
4939 !     real(kind=8) :: var,ene
4940
4941       open(icsa_seed,file=csa_seed,status="unknown")
4942        write(icsa_seed,*) "seed"
4943       close(31)
4944 #if defined(AIX) || defined(PGI)
4945        open(icsa_history,file=csa_history,status="unknown",&
4946         position="append")
4947 #else
4948        open(icsa_history,file=csa_history,status="unknown",&
4949         access="append")
4950 #endif
4951        write(icsa_history,*) nconf
4952        write(icsa_history,*) jstart,jend
4953        write(icsa_history,*) nstmax
4954        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4955        write(icsa_history,*) nran0,nran1,irr
4956        write(icsa_history,*) nseed
4957        write(icsa_history,*) ntotal,cut1,cut2
4958        write(icsa_history,*) estop
4959        write(icsa_history,*) icmax,irestart
4960        write(icsa_history,*) ntbankm,dele,difcut
4961        write(icsa_history,*) iref,rmscut,pnccut
4962        write(icsa_history,*) ndiff
4963
4964        write(icsa_history,*)
4965       close(icsa_history)
4966
4967       open(icsa_bank1,file=csa_bank1,status="unknown")
4968        write(icsa_bank1,*) 0
4969       close(icsa_bank1)
4970
4971       return
4972       end subroutine initial_write
4973 !-----------------------------------------------------------------------------
4974       subroutine restart_write
4975
4976       use csa_data
4977 !      implicit real*8 (a-h,o-z)
4978 !      include 'DIMENSIONS'
4979 !      include 'COMMON.IOUNITS'
4980 !      include 'COMMON.CSA'
4981 !      include 'COMMON.BANK'
4982 !el local variables
4983 !     integer :: ntf,ik,iw_pdb
4984 !     real(kind=8) :: var,ene
4985
4986 #if defined(AIX) || defined(PGI)
4987        open(icsa_history,file=csa_history,position="append")
4988 #else
4989        open(icsa_history,file=csa_history,access="append")
4990 #endif
4991        write(icsa_history,*)
4992        write(icsa_history,*) "This is restart"
4993        write(icsa_history,*)
4994        write(icsa_history,*) nconf
4995        write(icsa_history,*) jstart,jend
4996        write(icsa_history,*) nstmax
4997        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4998        write(icsa_history,*) nran0,nran1,irr
4999        write(icsa_history,*) nseed
5000        write(icsa_history,*) ntotal,cut1,cut2
5001        write(icsa_history,*) estop
5002        write(icsa_history,*) icmax,irestart
5003        write(icsa_history,*) ntbankm,dele,difcut
5004        write(icsa_history,*) iref,rmscut,pnccut
5005        write(icsa_history,*) ndiff
5006        write(icsa_history,*)
5007        write(icsa_history,*) "irestart is: ", irestart
5008
5009        write(icsa_history,*)
5010       close(icsa_history)
5011
5012       return
5013       end subroutine restart_write
5014 !-----------------------------------------------------------------------------
5015 ! test.F
5016 !-----------------------------------------------------------------------------
5017       subroutine write_pdb(npdb,titelloc,ee)
5018
5019 !      implicit real*8 (a-h,o-z)
5020 !      include 'DIMENSIONS'
5021 !      include 'COMMON.IOUNITS'
5022       character(len=50) :: titelloc1
5023       character*(*) :: titelloc
5024       character(len=3) :: zahl
5025       character(len=5) :: liczba5
5026       real(kind=8) :: ee
5027       integer :: npdb   !,ilen
5028 !el      external ilen
5029 !el local variables
5030       integer :: lenpre
5031 !     real(kind=8) :: var,ene
5032
5033       titelloc1=titelloc
5034       lenpre=ilen(prefix)
5035       if (npdb.lt.1000) then
5036        call numstr(npdb,zahl)
5037        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5038       else
5039         if (npdb.lt.10000) then
5040          write(liczba5,'(i1,i4)') 0,npdb
5041         else
5042          write(liczba5,'(i5)') npdb
5043         endif
5044         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5045       endif
5046       call pdbout(ee,titelloc1,ipdb)
5047       close(ipdb)
5048       return
5049       end subroutine write_pdb
5050 !-----------------------------------------------------------------------------
5051 ! thread.F
5052 !-----------------------------------------------------------------------------
5053       subroutine write_thread_summary
5054 ! Thread the sequence through a database of known structures
5055       use control_data, only: refstr
5056 !      use geometry
5057       use energy_data, only: n_ene_comp
5058       use compare_data
5059 !      implicit real*8 (a-h,o-z)
5060 !      include 'DIMENSIONS'
5061 #ifdef MPI
5062       use MPI_data      !include 'COMMON.INFO'
5063       include 'mpif.h'
5064 #endif
5065 !      include 'COMMON.CONTROL'
5066 !      include 'COMMON.CHAIN'
5067 !      include 'COMMON.DBASE'
5068 !      include 'COMMON.INTERACT'
5069 !      include 'COMMON.VAR'
5070 !      include 'COMMON.THREAD'
5071 !      include 'COMMON.FFIELD'
5072 !      include 'COMMON.SBRIDGE'
5073 !      include 'COMMON.HEADER'
5074 !      include 'COMMON.NAMES'
5075 !      include 'COMMON.IOUNITS'
5076 !      include 'COMMON.TIME1'
5077
5078       integer,dimension(maxthread) :: ip
5079       real(kind=8),dimension(0:n_ene) :: energia
5080 !el local variables
5081       integer :: i,j,ii,jj,ipj,ik,kk,ist
5082       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5083
5084       write (iout,'(30x,a/)') &
5085        '  *********** Summary threading statistics ************'
5086       write (iout,'(a)') 'Initial energies:'
5087       write (iout,'(a4,2x,a12,14a14,3a8)') &
5088        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5089        'RMSnat','NatCONT','NNCONT','RMS'
5090 ! Energy sort patterns
5091       do i=1,nthread
5092         ip(i)=i
5093       enddo
5094       do i=1,nthread-1
5095         enet=ener(n_ene-1,ip(i))
5096         jj=i
5097         do j=i+1,nthread
5098           if (ener(n_ene-1,ip(j)).lt.enet) then
5099             jj=j
5100             enet=ener(n_ene-1,ip(j))
5101           endif
5102         enddo
5103         if (jj.ne.i) then
5104           ipj=ip(jj)
5105           ip(jj)=ip(i)
5106           ip(i)=ipj
5107         endif
5108       enddo
5109       do ik=1,nthread
5110         i=ip(ik)
5111         ii=ipatt(1,i)
5112         ist=nres_base(2,ii)+ipatt(2,i)
5113         do kk=1,n_ene_comp
5114           energia(i)=ener0(kk,i)
5115         enddo
5116         etot=ener0(n_ene_comp+1,i)
5117         rmsnat=ener0(n_ene_comp+2,i)
5118         rms=ener0(n_ene_comp+3,i)
5119         frac=ener0(n_ene_comp+4,i)
5120         frac_nn=ener0(n_ene_comp+5,i)
5121
5122         if (refstr) then
5123         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5124         i,str_nam(ii),ist+1,&
5125         (energia(print_order(kk)),kk=1,nprint_ene),&
5126         etot,rmsnat,frac,frac_nn,rms
5127         else
5128         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5129         i,str_nam(ii),ist+1,&
5130         (energia(print_order(kk)),kk=1,nprint_ene),etot
5131         endif
5132       enddo
5133       write (iout,'(//a)') 'Final energies:'
5134       write (iout,'(a4,2x,a12,17a14,3a8)') &
5135        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5136        'RMSnat','NatCONT','NNCONT','RMS'
5137       do ik=1,nthread
5138         i=ip(ik)
5139         ii=ipatt(1,i)
5140         ist=nres_base(2,ii)+ipatt(2,i)
5141         do kk=1,n_ene_comp
5142           energia(kk)=ener(kk,ik)
5143         enddo
5144         etot=ener(n_ene_comp+1,i)
5145         rmsnat=ener(n_ene_comp+2,i)
5146         rms=ener(n_ene_comp+3,i)
5147         frac=ener(n_ene_comp+4,i)
5148         frac_nn=ener(n_ene_comp+5,i)
5149         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5150         i,str_nam(ii),ist+1,&
5151         (energia(print_order(kk)),kk=1,nprint_ene),&
5152         etot,rmsnat,frac,frac_nn,rms
5153       enddo
5154       write (iout,'(/a/)') 'IEXAM array:'
5155       write (iout,'(i5)') nexcl
5156       do i=1,nexcl
5157         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5158       enddo
5159       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5160        'Max. time for threading step ',max_time_for_thread,&
5161        'Average time for threading step: ',ave_time_for_thread
5162       return
5163       end subroutine write_thread_summary
5164 !-----------------------------------------------------------------------------
5165       subroutine write_stat_thread(ithread,ipattern,ist)
5166
5167       use energy_data, only: n_ene_comp
5168       use compare_data
5169 !      implicit real*8 (a-h,o-z)
5170 !      include "DIMENSIONS"
5171 !      include "COMMON.CONTROL"
5172 !      include "COMMON.IOUNITS"
5173 !      include "COMMON.THREAD"
5174 !      include "COMMON.FFIELD"
5175 !      include "COMMON.DBASE"
5176 !      include "COMMON.NAMES"
5177       real(kind=8),dimension(0:n_ene) :: energia
5178 !el local variables
5179       integer :: ithread,ipattern,ist,i
5180       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5181
5182 #if defined(AIX) || defined(PGI)
5183       open(istat,file=statname,position='append')
5184 #else
5185       open(istat,file=statname,access='append')
5186 #endif
5187       do i=1,n_ene_comp
5188         energia(i)=ener(i,ithread)
5189       enddo
5190       etot=ener(n_ene_comp+1,ithread)
5191       rmsnat=ener(n_ene_comp+2,ithread)
5192       rms=ener(n_ene_comp+3,ithread)
5193       frac=ener(n_ene_comp+4,ithread)
5194       frac_nn=ener(n_ene_comp+5,ithread)
5195       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5196         ithread,str_nam(ipattern),ist+1,&
5197         (energia(print_order(i)),i=1,nprint_ene),&
5198         etot,rmsnat,frac,frac_nn,rms
5199       close (istat)
5200       return
5201       end subroutine write_stat_thread
5202 !-----------------------------------------------------------------------------
5203 #endif
5204 !-----------------------------------------------------------------------------
5205       end module io_config