source/unres/io_config.f90

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