source/unres/io_config.f90

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