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