source/unres/io_config.f90

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