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