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