source/unres/io_config.f90

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