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)),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:maxtor,maxterm
 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
 747       integer :: ichir1,ichir2
 748 !      real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
 749 !el      allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
 750 !el      allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
 751 !
 752 ! For printing parameters after they are read set the following in the UNRES
 753 ! C-shell script:
 754 !
 755 ! setenv PRINT_PARM YES
 756 !
 757 ! To print parameters in LaTeX format rather than as ASCII tables:
 758 !
 759 ! setenv LATEX YES
 760 !
 761       call getenv_loc("PRINT_PARM",lancuch)
 762       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
 763       call getenv_loc("LATEX",lancuch)
 764       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
 765 !
 766       dwa16=2.0d0**(1.0d0/6.0d0)
 767       itypro=20
 768 ! Assign virtual-bond length
 769       vbl=3.8D0
 770       vblinv=1.0D0/vbl
 771       vblinv2=vblinv*vblinv
 772 !
 773 ! Read the virtual-bond parameters, masses, and moments of inertia
 774 ! and Stokes' radii of the peptide group and side chains
 775 !
 776       allocate(dsc(ntyp1)) !(ntyp1)
 777       allocate(dsc_inv(ntyp1)) !(ntyp1)
 778       allocate(nbondterm(ntyp)) !(ntyp)
 779       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 780       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 781       allocate(msc(ntyp+1)) !(ntyp+1)
 782       allocate(isc(ntyp+1)) !(ntyp+1)
 783       allocate(restok(ntyp+1)) !(ntyp+1)
 784       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
 785
 786       dsc(:)=0.0d0
 787       dsc_inv(:)=0.0d0
 788
 789 #ifdef CRYST_BOND
 790       read (ibond,*) vbldp0,akp,mp,ip,pstok
 791       do i=1,ntyp
 792         nbondterm(i)=1
 793         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
 794         dsc(i) = vbldsc0(1,i)
 795         if (i.eq.10) then
 796           dsc_inv(i)=0.0D0
 797         else
 798           dsc_inv(i)=1.0D0/dsc(i)
 799         endif
 800       enddo
 801 #else
 802       read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
 803       do i=1,ntyp
 804         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
 805          j=1,nbondterm(i)),msc(i),isc(i),restok(i)
 806         dsc(i) = vbldsc0(1,i)
 807         if (i.eq.10) then
 808           dsc_inv(i)=0.0D0
 809         else
 810           dsc_inv(i)=1.0D0/dsc(i)
 811         endif
 812       enddo
 813 #endif
 814       if (lprint) then
 815         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
 816         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
 817          'inertia','Pstok'
 818         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
 819         do i=1,ntyp
 820           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
 821             vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
 822           do j=2,nbondterm(i)
 823             write (iout,'(13x,3f10.5)') &
 824               vbldsc0(j,i),aksc(j,i),abond0(j,i)
 825           enddo
 826         enddo
 827       endif
 828 !----------------------------------------------------
 829       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
 830       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))      !(-ntyp:ntyp)
 831       allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
 832       allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
 833       allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
 834       allocate(gthet(3,-ntyp:ntyp))     !(3,-ntyp:ntyp)
 835       do i=-ntyp,ntyp
 836         a0thet(i)=0.0D0
 837         do j=1,2
 838          do ichir1=-1,1
 839           do ichir2=-1,1
 840           athet(j,i,ichir1,ichir2)=0.0D0
 841           bthet(j,i,ichir1,ichir2)=0.0D0
 842           enddo
 843          enddo
 844         enddo
 845         do j=0,3
 846           polthet(j,i)=0.0D0
 847         enddo
 848         do j=1,3
 849           gthet(j,i)=0.0D0
 850         enddo
 851         theta0(i)=0.0D0
 852         sig0(i)=0.0D0
 853         sigc0(i)=0.0D0
 854       enddo
 855
 856 #ifdef CRYST_THETA
 857 !
 858 ! Read the parameters of the probability distribution/energy expression
 859 ! of the virtual-bond valence angles theta
 860 !
 861       do i=1,ntyp
 862         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
 863           (bthet(j,i,1,1),j=1,2)
 864         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
 865         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
 866         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
 867         sigc0(i)=sigc0(i)**2
 868       enddo
 869       do i=1,ntyp
 870       athet(1,i,1,-1)=athet(1,i,1,1)
 871       athet(2,i,1,-1)=athet(2,i,1,1)
 872       bthet(1,i,1,-1)=-bthet(1,i,1,1)
 873       bthet(2,i,1,-1)=-bthet(2,i,1,1)
 874       athet(1,i,-1,1)=-athet(1,i,1,1)
 875       athet(2,i,-1,1)=-athet(2,i,1,1)
 876       bthet(1,i,-1,1)=bthet(1,i,1,1)
 877       bthet(2,i,-1,1)=bthet(2,i,1,1)
 878       enddo
 879       do i=-ntyp,-1
 880       a0thet(i)=a0thet(-i)
 881       athet(1,i,-1,-1)=athet(1,-i,1,1)
 882       athet(2,i,-1,-1)=-athet(2,-i,1,1)
 883       bthet(1,i,-1,-1)=bthet(1,-i,1,1)
 884       bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
 885       athet(1,i,-1,1)=athet(1,-i,1,1)
 886       athet(2,i,-1,1)=-athet(2,-i,1,1)
 887       bthet(1,i,-1,1)=-bthet(1,-i,1,1)
 888       bthet(2,i,-1,1)=bthet(2,-i,1,1)
 889       athet(1,i,1,-1)=-athet(1,-i,1,1)
 890       athet(2,i,1,-1)=athet(2,-i,1,1)
 891       bthet(1,i,1,-1)=bthet(1,-i,1,1)
 892       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
 893       theta0(i)=theta0(-i)
 894       sig0(i)=sig0(-i)
 895       sigc0(i)=sigc0(-i)
 896        do j=0,3
 897         polthet(j,i)=polthet(j,-i)
 898        enddo
 899        do j=1,3
 900          gthet(j,i)=gthet(j,-i)
 901        enddo
 902       enddo
 903
 904       close (ithep)
 905       if (lprint) then
 906       if (.not.LaTeX) then
 907         write (iout,'(a)') &
 908           'Parameters of the virtual-bond valence angles:'
 909         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
 910        '    ATHETA0   ','         A1   ','        A2    ',&
 911        '        B1    ','         B2   '
 912         do i=1,ntyp
 913           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
 914               a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
 915         enddo
 916         write (iout,'(/a/9x,5a/79(1h-))') &
 917        'Parameters of the expression for sigma(theta_c):',&
 918        '     ALPH0    ','      ALPH1   ','     ALPH2    ',&
 919        '     ALPH3    ','    SIGMA0C   '
 920         do i=1,ntyp
 921           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
 922             (polthet(j,i),j=0,3),sigc0(i)
 923         enddo
 924         write (iout,'(/a/9x,5a/79(1h-))') &
 925        'Parameters of the second gaussian:',&
 926        '    THETA0    ','     SIGMA0   ','        G1    ',&
 927        '        G2    ','         G3   '
 928         do i=1,ntyp
 929           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
 930              sig0(i),(gthet(j,i),j=1,3)
 931         enddo
 932        else
 933         write (iout,'(a)') &
 934           'Parameters of the virtual-bond valence angles:'
 935         write (iout,'(/a/9x,5a/79(1h-))') &
 936        'Coefficients of expansion',&
 937        '     theta0   ','    a1*10^2   ','   a2*10^2    ',&
 938        '   b1*10^1    ','    b2*10^1   '
 939         do i=1,ntyp
 940           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
 941               a0thet(i),(100*athet(j,i,1,1),j=1,2),&
 942               (10*bthet(j,i,1,1),j=1,2)
 943         enddo
 944         write (iout,'(/a/9x,5a/79(1h-))') &
 945        'Parameters of the expression for sigma(theta_c):',&
 946        ' alpha0       ','  alph1       ',' alph2        ',&
 947        ' alhp3        ','   sigma0c    '
 948         do i=1,ntyp
 949           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
 950             (polthet(j,i),j=0,3),sigc0(i)
 951         enddo
 952         write (iout,'(/a/9x,5a/79(1h-))') &
 953        'Parameters of the second gaussian:',&
 954        '    theta0    ','  sigma0*10^2 ','      G1*10^-1',&
 955        '        G2    ','   G3*10^1    '
 956         do i=1,ntyp
 957           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
 958              100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
 959         enddo
 960       endif
 961       endif
 962 #else
 963 !
 964 ! Read the parameters of Utheta determined from ab initio surfaces
 965 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
 966 !
 967       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
 968         ntheterm3,nsingle,ndouble
 969       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
 970
 971 !----------------------------------------------------
 972       allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
 973       allocate(aa0thet(-maxthetyp1:maxthetyp1,&
 974         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 975 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 976       allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,&
 977         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 978 !(maxtheterm,-maxthetyp1:maxthetyp1,&
 979 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 980       allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 981         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 982       allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 983         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 984       allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 985         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 986       allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
 987         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 988 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
 989 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
 990       allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
 991         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 992       allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
 993         -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 994 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
 995 !        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
 996
 997       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
 998       do i=-ntyp1,-1
 999         ithetyp(i)=-ithetyp(-i)
1000       enddo
1001       do iblock=1,2
1002       do i=-maxthetyp,maxthetyp
1003         do j=-maxthetyp,maxthetyp
1004           do k=-maxthetyp,maxthetyp
1005             aa0thet(i,j,k,iblock)=0.0d0
1006             do l=1,ntheterm
1007               aathet(l,i,j,k,iblock)=0.0d0
1008             enddo
1009             do l=1,ntheterm2
1010               do m=1,nsingle
1011                 bbthet(m,l,i,j,k,iblock)=0.0d0
1012                 ccthet(m,l,i,j,k,iblock)=0.0d0
1013                 ddthet(m,l,i,j,k,iblock)=0.0d0
1014                 eethet(m,l,i,j,k,iblock)=0.0d0
1015               enddo
1016             enddo
1017             do l=1,ntheterm3
1018               do m=1,ndouble
1019                 do mm=1,ndouble
1020                  ffthet(mm,m,l,i,j,k,iblock)=0.0d0
1021                  ggthet(mm,m,l,i,j,k,iblock)=0.0d0
1022                 enddo
1023               enddo
1024             enddo
1025           enddo
1026         enddo
1027       enddo
1028       enddo
1029 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1030       do iblock=1,2
1031 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1032 ! VAR:1=non-glicyne non-proline 2=proline
1033 ! VAR:negative values for D-aminoacid
1034       do i=0,nthetyp
1035         do j=-nthetyp,nthetyp
1036           do k=-nthetyp,nthetyp
1037             read (ithep,'(6a)',end=111,err=111) res1
1038             read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1039 ! VAR: aa0thet is variable describing the average value of Foureir
1040 ! VAR: expansion series
1041 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1042 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1043 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1044             read (ithep,*,end=111,err=111) &
1045               (aathet(l,i,j,k,iblock),l=1,ntheterm)
1046             read (ithep,*,end=111,err=111) &
1047              ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1048               (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1049               (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1050               (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1051               ll=1,ntheterm2)
1052             read (ithep,*,end=111,err=111) &
1053             (((ffthet(llll,lll,ll,i,j,k,iblock),&
1054                ffthet(lll,llll,ll,i,j,k,iblock),&
1055                ggthet(llll,lll,ll,i,j,k,iblock),&
1056                ggthet(lll,llll,ll,i,j,k,iblock),&
1057                llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1058           enddo
1059         enddo
1060       enddo
1061 !
1062 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1063 ! coefficients of theta-and-gamma-dependent terms are zero.
1064 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1065 ! RECOMENTDED AFTER VERSION 3.3)
1066 !      do i=1,nthetyp
1067 !        do j=1,nthetyp
1068 !          do l=1,ntheterm
1069 !            aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1070 !            aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1071 !          enddo
1072 !          aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1073 !          aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1074 !        enddo
1075 !        do l=1,ntheterm
1076 !          aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1077 !        enddo
1078 !        aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1079 !      enddo
1080 !      enddo
1081 ! AND COMMENT THE LOOPS BELOW
1082       do i=1,nthetyp
1083         do j=1,nthetyp
1084           do l=1,ntheterm
1085             aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1086             aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1087           enddo
1088           aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1089           aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1090         enddo
1091         do l=1,ntheterm
1092           aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1093         enddo
1094         aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1095       enddo
1096       enddo       !iblock
1097
1098 ! TILL HERE
1099 ! Substitution for D aminoacids from symmetry.
1100       do iblock=1,2
1101       do i=-nthetyp,0
1102         do j=-nthetyp,nthetyp
1103           do k=-nthetyp,nthetyp
1104            aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1105            do l=1,ntheterm
1106            aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1107            enddo
1108            do ll=1,ntheterm2
1109             do lll=1,nsingle
1110             bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1111             ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1112             ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1113             eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1114             enddo
1115           enddo
1116           do ll=1,ntheterm3
1117            do lll=2,ndouble
1118             do llll=1,lll-1
1119             ffthet(llll,lll,ll,i,j,k,iblock)= &
1120             ffthet(llll,lll,ll,-i,-j,-k,iblock)
1121             ffthet(lll,llll,ll,i,j,k,iblock)= &
1122             ffthet(lll,llll,ll,-i,-j,-k,iblock)
1123             ggthet(llll,lll,ll,i,j,k,iblock)= &
1124             -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1125             ggthet(lll,llll,ll,i,j,k,iblock)= &
1126             -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1127             enddo !ll
1128            enddo  !lll
1129           enddo   !llll
1130          enddo    !k
1131         enddo     !j
1132        enddo      !i
1133       enddo       !iblock
1134 !
1135 ! Control printout of the coefficients of virtual-bond-angle potentials
1136 !
1137       if (lprint) then
1138         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1139         do iblock=1,2
1140         do i=1,nthetyp+1
1141           do j=1,nthetyp+1
1142             do k=1,nthetyp+1
1143               write (iout,'(//4a)') &
1144                'Type ',onelett(i),onelett(j),onelett(k)
1145               write (iout,'(//a,10x,a)') " l","a[l]"
1146               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1147               write (iout,'(i2,1pe15.5)') &
1148                  (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1149             do l=1,ntheterm2
1150               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1151                 "b",l,"c",l,"d",l,"e",l
1152               do m=1,nsingle
1153                 write (iout,'(i2,4(1pe15.5))') m,&
1154                 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1155                 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1156               enddo
1157             enddo
1158             do l=1,ntheterm3
1159               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1160                 "f+",l,"f-",l,"g+",l,"g-",l
1161               do m=2,ndouble
1162                 do n=1,m-1
1163                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1164                     ffthet(n,m,l,i,j,k,iblock),&
1165                     ffthet(m,n,l,i,j,k,iblock),&
1166                     ggthet(n,m,l,i,j,k,iblock),&
1167                     ggthet(m,n,l,i,j,k,iblock)
1168                 enddo   !n
1169               enddo     !m
1170             enddo       !l
1171           enddo         !k
1172         enddo           !j
1173       enddo             !i
1174       enddo
1175       call flush(iout)
1176       endif
1177       write (2,*) "Start reading THETA_PDB",ithep_pdb
1178       do i=1,ntyp
1179 !      write (2,*) 'i=',i
1180         read (ithep_pdb,*,err=111,end=111) &
1181            a0thet(i),(athet(j,i,1,1),j=1,2),&
1182           (bthet(j,i,1,1),j=1,2)
1183         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1184         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1185         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1186         sigc0(i)=sigc0(i)**2
1187       enddo
1188       do i=1,ntyp
1189       athet(1,i,1,-1)=athet(1,i,1,1)
1190       athet(2,i,1,-1)=athet(2,i,1,1)
1191       bthet(1,i,1,-1)=-bthet(1,i,1,1)
1192       bthet(2,i,1,-1)=-bthet(2,i,1,1)
1193       athet(1,i,-1,1)=-athet(1,i,1,1)
1194       athet(2,i,-1,1)=-athet(2,i,1,1)
1195       bthet(1,i,-1,1)=bthet(1,i,1,1)
1196       bthet(2,i,-1,1)=bthet(2,i,1,1)
1197       enddo
1198       do i=-ntyp,-1
1199       a0thet(i)=a0thet(-i)
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       athet(1,i,1,-1)=-athet(1,-i,1,1)
1209       athet(2,i,1,-1)=athet(2,-i,1,1)
1210       bthet(1,i,1,-1)=bthet(1,-i,1,1)
1211       bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1212       theta0(i)=theta0(-i)
1213       sig0(i)=sig0(-i)
1214       sigc0(i)=sigc0(-i)
1215        do j=0,3
1216         polthet(j,i)=polthet(j,-i)
1217        enddo
1218        do j=1,3
1219          gthet(j,i)=gthet(j,-i)
1220        enddo
1221       enddo
1222       write (2,*) "End reading THETA_PDB"
1223       close (ithep_pdb)
1224 #endif
1225       close(ithep)
1226
1227 !-------------------------------------------
1228       allocate(nlob(ntyp1)) !(ntyp1)
1229       allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1230       allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1231       allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1232
1233       do i=1,ntyp
1234         do j=1,maxlob
1235           bsc(j,i)=0.0D0
1236           nlob(i)=0
1237         enddo
1238       enddo
1239       nlob(ntyp1)=0
1240       dsc(ntyp1)=0.0D0
1241
1242       do i=-ntyp,ntyp
1243         do j=1,maxlob
1244           do k=1,3
1245             censc(k,j,i)=0.0D0
1246           enddo
1247           do k=1,3
1248             do l=1,3
1249               gaussc(l,k,j,i)=0.0D0
1250             enddo
1251           enddo
1252         enddo
1253       enddo
1254
1255 #ifdef CRYST_SC
1256 !
1257 ! Read the parameters of the probability distribution/energy expression
1258 ! of the side chains.
1259 !
1260       do i=1,ntyp
1261         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1262         if (i.eq.10) then
1263           dsc_inv(i)=0.0D0
1264         else
1265           dsc_inv(i)=1.0D0/dsc(i)
1266         endif
1267         if (i.ne.10) then
1268         do j=1,nlob(i)
1269           do k=1,3
1270             do l=1,3
1271               blower(l,k,j)=0.0D0
1272             enddo
1273           enddo
1274         enddo
1275         bsc(1,i)=0.0D0
1276         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1277           ((blower(k,l,1),l=1,k),k=1,3)
1278         censc(1,1,-i)=censc(1,1,i)
1279         censc(2,1,-i)=censc(2,1,i)
1280         censc(3,1,-i)=-censc(3,1,i)
1281         do j=2,nlob(i)
1282           read (irotam,*,end=112,err=112) bsc(j,i)
1283           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1284                                        ((blower(k,l,j),l=1,k),k=1,3)
1285         censc(1,j,-i)=censc(1,j,i)
1286         censc(2,j,-i)=censc(2,j,i)
1287         censc(3,j,-i)=-censc(3,j,i)
1288 ! BSC is amplitude of Gaussian
1289         enddo
1290         do j=1,nlob(i)
1291           do k=1,3
1292             do l=1,k
1293               akl=0.0D0
1294               do m=1,3
1295                 akl=akl+blower(k,m,j)*blower(l,m,j)
1296               enddo
1297               gaussc(k,l,j,i)=akl
1298               gaussc(l,k,j,i)=akl
1299              if (((k.eq.3).and.(l.ne.3)) &
1300               .or.((l.eq.3).and.(k.ne.3))) then
1301                 gaussc(k,l,j,-i)=-akl
1302                 gaussc(l,k,j,-i)=-akl
1303               else
1304                 gaussc(k,l,j,-i)=akl
1305                 gaussc(l,k,j,-i)=akl
1306               endif
1307             enddo
1308           enddo
1309         enddo
1310         endif
1311       enddo
1312       close (irotam)
1313       if (lprint) then
1314         write (iout,'(/a)') 'Parameters of side-chain local geometry'
1315         do i=1,ntyp
1316           nlobi=nlob(i)
1317           if (nlobi.gt.0) then
1318             if (LaTeX) then
1319               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1320                ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1321                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1322                                    'log h',(bsc(j,i),j=1,nlobi)
1323                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1324               'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1325               do k=1,3
1326                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1327                        ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1328               enddo
1329             else
1330               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1331               write (iout,'(a,f10.4,4(16x,f10.4))') &
1332                                    'Center  ',(bsc(j,i),j=1,nlobi)
1333               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1334                  j=1,nlobi)
1335               write (iout,'(a)')
1336             endif
1337           endif
1338         enddo
1339       endif
1340 #else
1341 !
1342 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1343 ! added by Urszula Kozlowska 07/11/2007
1344 !
1345 !el Maximum number of SC local term fitting function coefficiants
1346 !el      integer,parameter :: maxsccoef=65
1347
1348       allocate(sc_parmin(65,ntyp))      !(maxsccoef,ntyp)
1349
1350       do i=1,ntyp
1351         read (irotam,*,end=112,err=112)
1352        if (i.eq.10) then
1353          read (irotam,*,end=112,err=112)
1354        else
1355          do j=1,65
1356            read(irotam,*,end=112,err=112) sc_parmin(j,i)
1357          enddo
1358        endif
1359       enddo
1360 !
1361 ! Read the parameters of the probability distribution/energy expression
1362 ! of the side chains.
1363 !
1364       write (2,*) "Start reading ROTAM_PDB"
1365       do i=1,ntyp
1366         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1367         if (i.eq.10) then
1368           dsc_inv(i)=0.0D0
1369         else
1370           dsc_inv(i)=1.0D0/dsc(i)
1371         endif
1372         if (i.ne.10) then
1373         do j=1,nlob(i)
1374           do k=1,3
1375             do l=1,3
1376               blower(l,k,j)=0.0D0
1377             enddo
1378           enddo
1379         enddo
1380         bsc(1,i)=0.0D0
1381         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1382           ((blower(k,l,1),l=1,k),k=1,3)
1383         do j=2,nlob(i)
1384           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1385           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1386                                        ((blower(k,l,j),l=1,k),k=1,3)
1387         enddo
1388         do j=1,nlob(i)
1389           do k=1,3
1390             do l=1,k
1391               akl=0.0D0
1392               do m=1,3
1393                 akl=akl+blower(k,m,j)*blower(l,m,j)
1394               enddo
1395               gaussc(k,l,j,i)=akl
1396               gaussc(l,k,j,i)=akl
1397             enddo
1398           enddo
1399         enddo
1400         endif
1401       enddo
1402       close (irotam_pdb)
1403       write (2,*) "End reading ROTAM_PDB"
1404 #endif
1405       close(irotam)
1406
1407 #ifdef CRYST_TOR
1408 !
1409 ! Read torsional parameters in old format
1410 !
1411       allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1412
1413       read (itorp,*,end=113,err=113) ntortyp,nterm_old
1414       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1415       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1416
1417 !el from energy module--------
1418       allocate(v1(nterm_old,ntortyp,ntortyp))
1419       allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1420 !el---------------------------
1421       do i=1,ntortyp
1422         do j=1,ntortyp
1423           read (itorp,'(a)')
1424           do k=1,nterm_old
1425             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1426           enddo
1427         enddo
1428       enddo
1429       close (itorp)
1430       if (lprint) then
1431         write (iout,'(/a/)') 'Torsional constants:'
1432         do i=1,ntortyp
1433           do j=1,ntortyp
1434             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1435             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1436           enddo
1437         enddo
1438       endif
1439 #else
1440 !
1441 ! Read torsional parameters
1442 !
1443       allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1444       read (itorp,*,end=113,err=113) ntortyp
1445 !el from energy module---------
1446       allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1447       allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1448
1449       allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1450       allocate(vlor2(maxlor,ntortyp,ntortyp))
1451       allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1452       allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1453
1454       allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1455       allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1456 !el---------------------------
1457       do iblock=1,2
1458         do i=-ntortyp,ntortyp
1459           do j=-ntortyp,ntortyp
1460             nterm(i,j,iblock)=0
1461             nlor(i,j,iblock)=0
1462           enddo
1463         enddo
1464       enddo
1465 !el---------------------------
1466
1467       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1468       do i=-ntyp,-1
1469        itortyp(i)=-itortyp(-i)
1470       enddo
1471 !      itortyp(ntyp1)=ntortyp
1472 !      itortyp(-ntyp1)=-ntortyp
1473       do iblock=1,2
1474       write (iout,*) 'ntortyp',ntortyp
1475       do i=0,ntortyp-1
1476         do j=-ntortyp+1,ntortyp-1
1477           read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1478                 nlor(i,j,iblock)
1479           nterm(-i,-j,iblock)=nterm(i,j,iblock)
1480           nlor(-i,-j,iblock)=nlor(i,j,iblock)
1481           v0ij=0.0d0
1482           si=-1.0d0
1483           do k=1,nterm(i,j,iblock)
1484             read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1485             v2(k,i,j,iblock)
1486             v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1487             v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1488             v0ij=v0ij+si*v1(k,i,j,iblock)
1489             si=-si
1490 !         write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1491 !         write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1492 !      v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1493           enddo
1494           do k=1,nlor(i,j,iblock)
1495             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1496               vlor2(k,i,j),vlor3(k,i,j)
1497             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1498           enddo
1499           v0(i,j,iblock)=v0ij
1500           v0(-i,-j,iblock)=v0ij
1501         enddo
1502       enddo
1503       enddo
1504       close (itorp)
1505       if (lprint) then
1506         write (iout,'(/a/)') 'Torsional constants:'
1507         do iblock=1,2
1508         do i=-ntortyp,ntortyp
1509           do j=-ntortyp,ntortyp
1510             write (iout,*) 'ityp',i,' jtyp',j
1511             write (iout,*) 'Fourier constants'
1512             do k=1,nterm(i,j,iblock)
1513               write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1514               v2(k,i,j,iblock)
1515             enddo
1516             write (iout,*) 'Lorenz constants'
1517             do k=1,nlor(i,j,iblock)
1518               write (iout,'(3(1pe15.5))') &
1519                vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1520             enddo
1521           enddo
1522         enddo
1523         enddo
1524       endif
1525 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1526 !
1527 ! 6/23/01 Read parameters for double torsionals
1528 !
1529 !el from energy module------------
1530       allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1531       allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1532 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1533       allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1534       allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1535         !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1536       allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1537       allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1538         !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1539 !---------------------------------
1540
1541       do iblock=1,2
1542       do i=0,ntortyp-1
1543         do j=-ntortyp+1,ntortyp-1
1544           do k=-ntortyp+1,ntortyp-1
1545             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1546 !              write (iout,*) "OK onelett",
1547 !     &         i,j,k,t1,t2,t3
1548
1549             if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1550               .or. t3.ne.toronelet(k)) then
1551              write (iout,*) "Error in double torsional parameter file",&
1552                i,j,k,t1,t2,t3
1553 #ifdef MPI
1554               call MPI_Finalize(Ierror)
1555 #endif
1556                stop "Error in double torsional parameter file"
1557             endif
1558            read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1559                ntermd_2(i,j,k,iblock)
1560             ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1561             ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1562             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1563                ntermd_1(i,j,k,iblock))
1564             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1565                ntermd_1(i,j,k,iblock))
1566             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1567                ntermd_1(i,j,k,iblock))
1568             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1569                ntermd_1(i,j,k,iblock))
1570 ! Martix of D parameters for one dimesional foureir series
1571             do l=1,ntermd_1(i,j,k,iblock)
1572              v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1573              v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1574              v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1575              v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1576 !            write(iout,*) "whcodze" ,
1577 !     & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1578             enddo
1579             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1580                v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1581                v2s(m,l,i,j,k,iblock),&
1582                m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1583 ! Martix of D parameters for two dimesional fourier series
1584             do l=1,ntermd_2(i,j,k,iblock)
1585              do m=1,l-1
1586              v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1587              v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1588              v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1589              v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1590              enddo!m
1591             enddo!l
1592           enddo!k
1593         enddo!j
1594       enddo!i
1595       enddo!iblock
1596       if (lprint) then
1597       write (iout,*)
1598       write (iout,*) 'Constants for double torsionals'
1599       do iblock=1,2
1600       do i=0,ntortyp-1
1601         do j=-ntortyp+1,ntortyp-1
1602           do k=-ntortyp+1,ntortyp-1
1603             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1604               ' nsingle',ntermd_1(i,j,k,iblock),&
1605               ' ndouble',ntermd_2(i,j,k,iblock)
1606             write (iout,*)
1607             write (iout,*) 'Single angles:'
1608             do l=1,ntermd_1(i,j,k,iblock)
1609               write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1610                  v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1611                  v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1612                  v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1613             enddo
1614             write (iout,*)
1615             write (iout,*) 'Pairs of angles:'
1616             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1617             do l=1,ntermd_2(i,j,k,iblock)
1618               write (iout,'(i5,20f10.5)') &
1619                l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1620             enddo
1621             write (iout,*)
1622             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1623             do l=1,ntermd_2(i,j,k,iblock)
1624               write (iout,'(i5,20f10.5)') &
1625                l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1626                (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1627             enddo
1628             write (iout,*)
1629           enddo
1630         enddo
1631       enddo
1632       enddo
1633       endif
1634 #endif
1635 ! Read of Side-chain backbone correlation parameters
1636 ! Modified 11 May 2012 by Adasko
1637 !CC
1638 !
1639       read (isccor,*,end=119,err=119) nsccortyp
1640
1641 !el from module energy-------------
1642       allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1643       allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1644       allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1645       allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp))   !(maxterm_sccor,20,20)
1646 !-----------------------------------
1647 #ifdef SCCORPDB
1648 !el from module energy-------------
1649       allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1650
1651       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1652       do i=-ntyp,-1
1653         isccortyp(i)=-isccortyp(-i)
1654       enddo
1655       iscprol=isccortyp(20)
1656 !      write (iout,*) 'ntortyp',ntortyp
1657       maxinter=3
1658 !c maxinter is maximum interaction sites
1659 !el from module energy---------
1660       allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1661       allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1662                -nsccortyp:nsccortyp))
1663       allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1664                -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1665       allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1666                -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1667 !-----------------------------------
1668       do l=1,maxinter
1669       do i=1,nsccortyp
1670         do j=1,nsccortyp
1671           read (isccor,*,end=119,err=119) &
1672       nterm_sccor(i,j),nlor_sccor(i,j)
1673           v0ijsccor=0.0d0
1674           v0ijsccor1=0.0d0
1675           v0ijsccor2=0.0d0
1676           v0ijsccor3=0.0d0
1677           si=-1.0d0
1678           nterm_sccor(-i,j)=nterm_sccor(i,j)
1679           nterm_sccor(-i,-j)=nterm_sccor(i,j)
1680           nterm_sccor(i,-j)=nterm_sccor(i,j)
1681           do k=1,nterm_sccor(i,j)
1682             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1683            v2sccor(k,l,i,j)
1684             if (j.eq.iscprol) then
1685              if (i.eq.isccortyp(10)) then
1686              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1687              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1688              else
1689              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1690                               +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1691              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1692                               +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1693              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1694              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1695              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1696              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1697              endif
1698             else
1699              if (i.eq.isccortyp(10)) then
1700              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1701              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1702              else
1703                if (j.eq.isccortyp(10)) then
1704              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1705              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1706                else
1707              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1708              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1709              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1710              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1711              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1712              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1713                 endif
1714                endif
1715             endif
1716             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1717             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1718             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1719             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1720             si=-si
1721           enddo
1722           do k=1,nlor_sccor(i,j)
1723             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1724               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1725             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1726       (1+vlor3sccor(k,i,j)**2)
1727           enddo
1728           v0sccor(l,i,j)=v0ijsccor
1729           v0sccor(l,-i,j)=v0ijsccor1
1730           v0sccor(l,i,-j)=v0ijsccor2
1731           v0sccor(l,-i,-j)=v0ijsccor3
1732         enddo
1733       enddo
1734       enddo
1735       close (isccor)
1736 #else
1737 !el from module energy-------------
1738       allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1739
1740       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1741 !      write (iout,*) 'ntortyp',ntortyp
1742       maxinter=3
1743 !c maxinter is maximum interaction sites
1744 !el from module energy---------
1745       allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1746       allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1747       allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1748       allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1749 !-----------------------------------
1750       do l=1,maxinter
1751       do i=1,nsccortyp
1752         do j=1,nsccortyp
1753           read (isccor,*,end=119,err=119) &
1754        nterm_sccor(i,j),nlor_sccor(i,j)
1755           v0ijsccor=0.0d0
1756           si=-1.0d0
1757
1758           do k=1,nterm_sccor(i,j)
1759             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1760            v2sccor(k,l,i,j)
1761             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1762             si=-si
1763           enddo
1764           do k=1,nlor_sccor(i,j)
1765             read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1766               vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1767             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1768       (1+vlor3sccor(k,i,j)**2)
1769           enddo
1770           v0sccor(l,i,j)=v0ijsccor !el ,iblock
1771         enddo
1772       enddo
1773       enddo
1774       close (isccor)
1775
1776 #endif
1777       if (lprint) then
1778         write (iout,'(/a/)') 'Torsional constants:'
1779         do i=1,nsccortyp
1780           do j=1,nsccortyp
1781             write (iout,*) 'ityp',i,' jtyp',j
1782             write (iout,*) 'Fourier constants'
1783             do k=1,nterm_sccor(i,j)
1784       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1785             enddo
1786             write (iout,*) 'Lorenz constants'
1787             do k=1,nlor_sccor(i,j)
1788               write (iout,'(3(1pe15.5))') &
1789                vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1790             enddo
1791           enddo
1792         enddo
1793       endif
1794
1795 !
1796 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1797 !         interaction energy of the Gly, Ala, and Pro prototypes.
1798 !
1799
1800       if (lprint) then
1801         write (iout,*)
1802         write (iout,*) "Coefficients of the cumulants"
1803       endif
1804       read (ifourier,*) nloctyp
1805 !write(iout,*) "nloctyp",nloctyp
1806 !el from module energy-------
1807       allocate(b1(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1808       allocate(b2(2,-nloctyp-1:nloctyp+1))      !(2,-maxtor:maxtor)
1809       allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1810       allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1811       allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1812       allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1813       allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1814       allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1815 ! el
1816         b1(1,:)=0.0d0
1817         b1(2,:)=0.0d0
1818 !--------------------------------
1819
1820       do i=0,nloctyp-1
1821         read (ifourier,*,end=115,err=115)
1822         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1823         if (lprint) then
1824           write (iout,*) 'Type',i
1825           write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1826         endif
1827         B1(1,i)  = b(3)
1828         B1(2,i)  = b(5)
1829         B1(1,-i) = b(3)
1830         B1(2,-i) = -b(5)
1831 !        b1(1,i)=0.0d0
1832 !        b1(2,i)=0.0d0
1833         B1tilde(1,i) = b(3)
1834         B1tilde(2,i) =-b(5)
1835         B1tilde(1,-i) =-b(3)
1836         B1tilde(2,-i) =b(5)
1837 !        b1tilde(1,i)=0.0d0
1838 !        b1tilde(2,i)=0.0d0
1839         B2(1,i)  = b(2)
1840         B2(2,i)  = b(4)
1841         B2(1,-i)  =b(2)
1842         B2(2,-i)  =-b(4)
1843
1844 !        b2(1,i)=0.0d0
1845 !        b2(2,i)=0.0d0
1846         CC(1,1,i)= b(7)
1847         CC(2,2,i)=-b(7)
1848         CC(2,1,i)= b(9)
1849         CC(1,2,i)= b(9)
1850         CC(1,1,-i)= b(7)
1851         CC(2,2,-i)=-b(7)
1852         CC(2,1,-i)=-b(9)
1853         CC(1,2,-i)=-b(9)
1854 !        CC(1,1,i)=0.0d0
1855 !        CC(2,2,i)=0.0d0
1856 !        CC(2,1,i)=0.0d0
1857 !        CC(1,2,i)=0.0d0
1858         Ctilde(1,1,i)=b(7)
1859         Ctilde(1,2,i)=b(9)
1860         Ctilde(2,1,i)=-b(9)
1861         Ctilde(2,2,i)=b(7)
1862         Ctilde(1,1,-i)=b(7)
1863         Ctilde(1,2,-i)=-b(9)
1864         Ctilde(2,1,-i)=b(9)
1865         Ctilde(2,2,-i)=b(7)
1866
1867 !        Ctilde(1,1,i)=0.0d0
1868 !        Ctilde(1,2,i)=0.0d0
1869 !        Ctilde(2,1,i)=0.0d0
1870 !        Ctilde(2,2,i)=0.0d0
1871         DD(1,1,i)= b(6)
1872         DD(2,2,i)=-b(6)
1873         DD(2,1,i)= b(8)
1874         DD(1,2,i)= b(8)
1875         DD(1,1,-i)= b(6)
1876         DD(2,2,-i)=-b(6)
1877         DD(2,1,-i)=-b(8)
1878         DD(1,2,-i)=-b(8)
1879 !        DD(1,1,i)=0.0d0
1880 !        DD(2,2,i)=0.0d0
1881 !        DD(2,1,i)=0.0d0
1882 !        DD(1,2,i)=0.0d0
1883         Dtilde(1,1,i)=b(6)
1884         Dtilde(1,2,i)=b(8)
1885         Dtilde(2,1,i)=-b(8)
1886         Dtilde(2,2,i)=b(6)
1887         Dtilde(1,1,-i)=b(6)
1888         Dtilde(1,2,-i)=-b(8)
1889         Dtilde(2,1,-i)=b(8)
1890         Dtilde(2,2,-i)=b(6)
1891
1892 !        Dtilde(1,1,i)=0.0d0
1893 !        Dtilde(1,2,i)=0.0d0
1894 !        Dtilde(2,1,i)=0.0d0
1895 !        Dtilde(2,2,i)=0.0d0
1896         EE(1,1,i)= b(10)+b(11)
1897         EE(2,2,i)=-b(10)+b(11)
1898         EE(2,1,i)= b(12)-b(13)
1899         EE(1,2,i)= b(12)+b(13)
1900         EE(1,1,-i)= b(10)+b(11)
1901         EE(2,2,-i)=-b(10)+b(11)
1902         EE(2,1,-i)=-b(12)+b(13)
1903         EE(1,2,-i)=-b(12)-b(13)
1904
1905 !        ee(1,1,i)=1.0d0
1906 !        ee(2,2,i)=1.0d0
1907 !        ee(2,1,i)=0.0d0
1908 !        ee(1,2,i)=0.0d0
1909 !        ee(2,1,i)=ee(1,2,i)
1910       enddo
1911       if (lprint) then
1912       do i=1,nloctyp
1913         write (iout,*) 'Type',i
1914         write (iout,*) 'B1'
1915         write(iout,*) B1(1,i),B1(2,i)
1916         write (iout,*) 'B2'
1917         write(iout,*) B2(1,i),B2(2,i)
1918         write (iout,*) 'CC'
1919         do j=1,2
1920           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1921         enddo
1922         write(iout,*) 'DD'
1923         do j=1,2
1924           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1925         enddo
1926         write(iout,*) 'EE'
1927         do j=1,2
1928           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1929         enddo
1930       enddo
1931       endif
1932 !
1933 ! Read electrostatic-interaction parameters
1934 !
1935
1936       if (lprint) then
1937         write (iout,*)
1938         write (iout,'(/a)') 'Electrostatic interaction constants:'
1939         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1940                   'IT','JT','APP','BPP','AEL6','AEL3'
1941       endif
1942       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1943       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1944       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1945       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1946       close (ielep)
1947       do i=1,2
1948         do j=1,2
1949         rri=rpp(i,j)**6
1950         app (i,j)=epp(i,j)*rri*rri
1951         bpp (i,j)=-2.0D0*epp(i,j)*rri
1952         ael6(i,j)=elpp6(i,j)*4.2D0**6
1953         ael3(i,j)=elpp3(i,j)*4.2D0**3
1954 !        lprint=.true.
1955         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1956                           ael6(i,j),ael3(i,j)
1957 !        lprint=.false.
1958         enddo
1959       enddo
1960 !
1961 ! Read side-chain interaction parameters.
1962 !
1963 !el from module energy - COMMON.INTERACT-------
1964       allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1965       allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1966       allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1967       allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1968       allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1969       do i=1,ntyp
1970         do j=1,ntyp
1971           augm(i,j)=0.0D0
1972         enddo
1973         chip(i)=0.0D0
1974         alp(i)=0.0D0
1975         sigma0(i)=0.0D0
1976         sigii(i)=0.0D0
1977         rr0(i)=0.0D0
1978       enddo
1979 !--------------------------------
1980
1981       read (isidep,*,end=117,err=117) ipot,expon
1982       if (ipot.lt.1 .or. ipot.gt.5) then
1983         write (iout,'(2a)') 'Error while reading SC interaction',&
1984                      'potential file - unknown potential type.'
1985 #ifdef MPI
1986         call MPI_Finalize(Ierror)
1987 #endif
1988         stop
1989       endif
1990       expon2=expon/2
1991       if(me.eq.king) &
1992        write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1993        ', exponents are ',expon,2*expon
1994 !      goto (10,20,30,30,40) ipot
1995       select case(ipot)
1996 !----------------------- LJ potential ---------------------------------
1997        case (1)
1998 !   10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1999          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2000            (sigma0(i),i=1,ntyp)
2001         if (lprint) then
2002           write (iout,'(/a/)') 'Parameters of the LJ potential:'
2003           write (iout,'(a/)') 'The epsilon array:'
2004           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2005           write (iout,'(/a)') 'One-body parameters:'
2006           write (iout,'(a,4x,a)') 'residue','sigma'
2007           write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
2008         endif
2009 !      goto 50
2010 !----------------------- LJK potential --------------------------------
2011        case(2)
2012 !   20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2013          read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2014           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2015         if (lprint) then
2016           write (iout,'(/a/)') 'Parameters of the LJK potential:'
2017           write (iout,'(a/)') 'The epsilon array:'
2018           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2019           write (iout,'(/a)') 'One-body parameters:'
2020           write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
2021           write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
2022                 i=1,ntyp)
2023         endif
2024 !      goto 50
2025 !---------------------- GB or BP potential -----------------------------
2026        case(3:4)
2027 !   30 do i=1,ntyp
2028         do i=1,ntyp
2029          read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2030         enddo
2031         read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2032         read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2033         read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2034         read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2035 ! For the GB potential convert sigma'**2 into chi'
2036         if (ipot.eq.4) then
2037           do i=1,ntyp
2038             chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2039           enddo
2040         endif
2041         if (lprint) then
2042           write (iout,'(/a/)') 'Parameters of the BP potential:'
2043           write (iout,'(a/)') 'The epsilon array:'
2044           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2045           write (iout,'(/a)') 'One-body parameters:'
2046           write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',&
2047                '    chip  ','    alph  '
2048           write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
2049                              chip(i),alp(i),i=1,ntyp)
2050         endif
2051 !      goto 50
2052 !--------------------- GBV potential -----------------------------------
2053        case(5)
2054 !   40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2055         read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2056           (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2057           (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2058         if (lprint) then
2059           write (iout,'(/a/)') 'Parameters of the GBV potential:'
2060           write (iout,'(a/)') 'The epsilon array:'
2061           call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2062           write (iout,'(/a)') 'One-body parameters:'
2063           write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',&
2064               's||/s_|_^2','    chip  ','    alph  '
2065           write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2066                    sigii(i),chip(i),alp(i),i=1,ntyp)
2067         endif
2068        case default
2069         write(iout,*)"Wrong ipot"
2070 !   50 continue
2071       end select
2072       continue
2073       close (isidep)
2074 !-----------------------------------------------------------------------
2075 ! Calculate the "working" parameters of SC interactions.
2076
2077 !el from module energy - COMMON.INTERACT-------
2078       allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2079       allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2080       do i=1,ntyp1
2081         do j=1,ntyp1
2082           aa(i,j)=0.0D0
2083           bb(i,j)=0.0D0
2084           chi(i,j)=0.0D0
2085           sigma(i,j)=0.0D0
2086           r0(i,j)=0.0D0
2087         enddo
2088       enddo
2089 !--------------------------------
2090
2091       do i=2,ntyp
2092         do j=1,i-1
2093           eps(i,j)=eps(j,i)
2094         enddo
2095       enddo
2096       do i=1,ntyp
2097         do j=i,ntyp
2098           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2099           sigma(j,i)=sigma(i,j)
2100           rs0(i,j)=dwa16*sigma(i,j)
2101           rs0(j,i)=rs0(i,j)
2102         enddo
2103       enddo
2104       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2105        'Working parameters of the SC interactions:',&
2106        '     a    ','     b    ','   augm   ','  sigma ','   r0   ',&
2107        '  chi1   ','   chi2   '
2108       do i=1,ntyp
2109         do j=i,ntyp
2110           epsij=eps(i,j)
2111           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2112             rrij=sigma(i,j)
2113           else
2114             rrij=rr0(i)+rr0(j)
2115           endif
2116           r0(i,j)=rrij
2117           r0(j,i)=rrij
2118           rrij=rrij**expon
2119           epsij=eps(i,j)
2120           sigeps=dsign(1.0D0,epsij)
2121           epsij=dabs(epsij)
2122           aa(i,j)=epsij*rrij*rrij
2123           bb(i,j)=-sigeps*epsij*rrij
2124           aa(j,i)=aa(i,j)
2125           bb(j,i)=bb(i,j)
2126           if (ipot.gt.2) then
2127             sigt1sq=sigma0(i)**2
2128             sigt2sq=sigma0(j)**2
2129             sigii1=sigii(i)
2130             sigii2=sigii(j)
2131             ratsig1=sigt2sq/sigt1sq
2132             ratsig2=1.0D0/ratsig1
2133             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2134             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2135             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2136           else
2137             rsum_max=sigma(i,j)
2138           endif
2139 !         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2140             sigmaii(i,j)=rsum_max
2141             sigmaii(j,i)=rsum_max
2142 !         else
2143 !           sigmaii(i,j)=r0(i,j)
2144 !           sigmaii(j,i)=r0(i,j)
2145 !         endif
2146 !d        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2147           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2148             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2149             augm(i,j)=epsij*r_augm**(2*expon)
2150 !           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2151             augm(j,i)=augm(i,j)
2152           else
2153             augm(i,j)=0.0D0
2154             augm(j,i)=0.0D0
2155           endif
2156           if (lprint) then
2157             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2158             restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
2159             sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2160           endif
2161         enddo
2162       enddo
2163
2164       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2165       do i=1,ntyp
2166         do j=1,2
2167           bad(i,j)=0.0D0
2168         enddo
2169       enddo
2170
2171 #ifdef OLDSCP
2172 !
2173 ! Define the SC-p interaction constants (hard-coded; old style)
2174 !
2175       do i=1,ntyp
2176 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2177 ! helix formation)
2178 !       aad(i,1)=0.3D0*4.0D0**12
2179 ! Following line for constants currently implemented
2180 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2181         aad(i,1)=1.5D0*4.0D0**12
2182 !       aad(i,1)=0.17D0*5.6D0**12
2183         aad(i,2)=aad(i,1)
2184 ! "Soft" SC-p repulsion
2185         bad(i,1)=0.0D0
2186 ! Following line for constants currently implemented
2187 !       aad(i,1)=0.3D0*4.0D0**6
2188 ! "Hard" SC-p repulsion
2189         bad(i,1)=3.0D0*4.0D0**6
2190 !       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2191         bad(i,2)=bad(i,1)
2192 !       aad(i,1)=0.0D0
2193 !       aad(i,2)=0.0D0
2194 !       bad(i,1)=1228.8D0
2195 !       bad(i,2)=1228.8D0
2196       enddo
2197 #else
2198 !
2199 ! 8/9/01 Read the SC-p interaction constants from file
2200 !
2201       do i=1,ntyp
2202         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2203       enddo
2204       do i=1,ntyp
2205         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2206         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2207         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2208         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2209       enddo
2210 !      lprint=.true.
2211       if (lprint) then
2212         write (iout,*) "Parameters of SC-p interactions:"
2213         do i=1,ntyp
2214           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2215            eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2216         enddo
2217       endif
2218 !      lprint=.false.
2219 #endif
2220 !
2221 ! Define the constants of the disulfide bridge
2222 !
2223       ebr=-5.50D0
2224 !
2225 ! Old arbitrary potential - commented out.
2226 !
2227 !      dbr= 4.20D0
2228 !      fbr= 3.30D0
2229 !
2230 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2231 ! energy surface of diethyl disulfide.
2232 ! A. Liwo and U. Kozlowska, 11/24/03
2233 !
2234       D0CM = 3.78d0
2235       AKCM = 15.1d0
2236       AKTH = 11.0d0
2237       AKCT = 12.0d0
2238       V1SS =-1.08d0
2239       V2SS = 7.61d0
2240       V3SS = 13.7d0
2241 !      akcm=0.0d0
2242 !      akth=0.0d0
2243 !      akct=0.0d0
2244 !      v1ss=0.0d0
2245 !      v2ss=0.0d0
2246 !      v3ss=0.0d0
2247
2248       if(me.eq.king) then
2249       write (iout,'(/a)') "Disulfide bridge parameters:"
2250       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2251       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2252       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2253       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2254         ' v3ss:',v3ss
2255       endif
2256       return
2257   111 write (iout,*) "Error reading bending energy parameters."
2258       goto 999
2259   112 write (iout,*) "Error reading rotamer energy parameters."
2260       goto 999
2261   113 write (iout,*) "Error reading torsional energy parameters."
2262       goto 999
2263   114 write (iout,*) "Error reading double torsional energy parameters."
2264       goto 999
2265   115 write (iout,*) &
2266         "Error reading cumulant (multibody energy) parameters."
2267       goto 999
2268   116 write (iout,*) "Error reading electrostatic energy parameters."
2269       goto 999
2270   117 write (iout,*) "Error reading side chain interaction parameters."
2271       goto 999
2272   118 write (iout,*) "Error reading SCp interaction parameters."
2273       goto 999
2274   119 write (iout,*) "Error reading SCCOR parameters"
2275   999 continue
2276 #ifdef MPI
2277       call MPI_Finalize(Ierror)
2278 #endif
2279       stop
2280       return
2281       end subroutine parmread
2282 #endif
2283 !-----------------------------------------------------------------------------
2284 ! printmat.f
2285 !-----------------------------------------------------------------------------
2286       subroutine printmat(ldim,m,n,iout,key,a)
2287
2288       integer :: n,ldim
2289       character(len=3),dimension(n) :: key
2290       real(kind=8),dimension(ldim,n) :: a
2291 !el local variables
2292       integer :: i,j,k,m,iout,nlim
2293
2294       do 1 i=1,n,8
2295       nlim=min0(i+7,n)
2296       write (iout,1000) (key(k),k=i,nlim)
2297       write (iout,1020)
2298  1000 format (/5x,8(6x,a3))
2299  1020 format (/80(1h-)/)
2300       do 2 j=1,n
2301       write (iout,1010) key(j),(a(j,k),k=i,nlim)
2302     2 continue
2303     1 continue
2304  1010 format (a3,2x,8(f9.4))
2305       return
2306       end subroutine printmat
2307 !-----------------------------------------------------------------------------
2308 ! readpdb.F
2309 !-----------------------------------------------------------------------------
2310       subroutine readpdb
2311 ! Read the PDB file and convert the peptide geometry into virtual-chain
2312 ! geometry.
2313       use geometry_data
2314       use energy_data, only: itype
2315       use control_data
2316       use compare_data
2317       use MPI_data
2318       use control, only: rescode
2319 !      implicit real*8 (a-h,o-z)
2320 !      include 'DIMENSIONS'
2321 !      include 'COMMON.LOCAL'
2322 !      include 'COMMON.VAR'
2323 !      include 'COMMON.CHAIN'
2324 !      include 'COMMON.INTERACT'
2325 !      include 'COMMON.IOUNITS'
2326 !      include 'COMMON.GEO'
2327 !      include 'COMMON.NAMES'
2328 !      include 'COMMON.CONTROL'
2329 !      include 'COMMON.DISTFIT'
2330 !      include 'COMMON.SETUP'
2331       integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2332 !        ishift_pdb
2333       logical :: lprn=.true.,fail
2334       real(kind=8),dimension(3) :: e1,e2,e3
2335       real(kind=8) :: dcj,efree_temp
2336       character(len=3) :: seq,res
2337       character(len=5) :: atom
2338       character(len=80) :: card
2339       real(kind=8),dimension(3,20) :: sccor
2340       integer :: kkk,lll,icha,kupa      !rescode,
2341       real(kind=8) :: cou
2342 !el local varables
2343       integer,dimension(2,maxres/3) :: hfrag_alloc
2344       integer,dimension(4,maxres/3) :: bfrag_alloc
2345       real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2346
2347       efree_temp=0.0d0
2348       ibeg=1
2349       ishift1=0
2350       ishift=0
2351 !      write (2,*) "UNRES_PDB",unres_pdb
2352       ires=0
2353       ires_old=0
2354       nres=0
2355       iii=0
2356       lsecondary=.false.
2357       nhfrag=0
2358       nbfrag=0
2359 !-----------------------------
2360       allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2361       allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2362
2363       do i=1,100000
2364         read (ipdbin,'(a80)',end=10) card
2365 !       write (iout,'(a)') card
2366         if (card(:5).eq.'HELIX') then
2367           nhfrag=nhfrag+1
2368           lsecondary=.true.
2369           read(card(22:25),*) hfrag(1,nhfrag)
2370           read(card(34:37),*) hfrag(2,nhfrag)
2371         endif
2372         if (card(:5).eq.'SHEET') then
2373           nbfrag=nbfrag+1
2374           lsecondary=.true.
2375           read(card(24:26),*) bfrag(1,nbfrag)
2376           read(card(35:37),*) bfrag(2,nbfrag)
2377 !rc----------------------------------------
2378 !rc  to be corrected !!!
2379           bfrag(3,nbfrag)=bfrag(1,nbfrag)
2380           bfrag(4,nbfrag)=bfrag(2,nbfrag)
2381 !rc----------------------------------------
2382         endif
2383         if (card(:3).eq.'END') then
2384           goto 10
2385         else if (card(:3).eq.'TER') then
2386 ! End current chain
2387           ires_old=ires+1
2388           ishift1=ishift1+1
2389           itype(ires_old)=ntyp1
2390           ibeg=2
2391 !          write (iout,*) "Chain ended",ires,ishift,ires_old
2392           if (unres_pdb) then
2393             do j=1,3
2394               dc(j,ires)=sccor(j,iii)
2395             enddo
2396           else
2397             call sccenter(ires,iii,sccor)
2398           endif
2399           iii=0
2400         endif
2401 ! Read free energy
2402         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2403 ! Fish out the ATOM cards.
2404         if (index(card(1:4),'ATOM').gt.0) then
2405           read (card(12:16),*) atom
2406 !          write (iout,*) "! ",atom," !",ires
2407 !          if (atom.eq.'CA' .or. atom.eq.'CH3') then
2408           read (card(23:26),*) ires
2409           read (card(18:20),'(a3)') res
2410 !          write (iout,*) "ires",ires,ires-ishift+ishift1,
2411 !     &      " ires_old",ires_old
2412 !          write (iout,*) "ishift",ishift," ishift1",ishift1
2413 !          write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2414           if (ires-ishift+ishift1.ne.ires_old) then
2415 ! Calculate the CM of the preceding residue.
2416 !            if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2417             if (ibeg.eq.0) then
2418 !              write (iout,*) "Calculating sidechain center iii",iii
2419               if (unres_pdb) then
2420                 do j=1,3
2421                   dc(j,ires+nres)=sccor(j,iii)
2422                 enddo
2423               else
2424                 call sccenter(ires_old,iii,sccor)
2425               endif
2426               iii=0
2427             endif
2428 ! Start new residue.
2429             if (res.eq.'Cl-' .or. res.eq.'Na+') then
2430               ires=ires_old
2431               cycle
2432             else if (ibeg.eq.1) then
2433               write (iout,*) "BEG ires",ires
2434               ishift=ires-1
2435               if (res.ne.'GLY' .and. res.ne. 'ACE') then
2436                 ishift=ishift-1
2437                 itype(1)=ntyp1
2438               endif
2439               ires=ires-ishift+ishift1
2440               ires_old=ires
2441 !              write (iout,*) "ishift",ishift," ires",ires,&
2442 !               " ires_old",ires_old
2443               ibeg=0
2444             else if (ibeg.eq.2) then
2445 ! Start a new chain
2446               ishift=-ires_old+ires-1 !!!!!
2447               ishift1=ishift1-1    !!!!!
2448 !              write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2449               ires=ires-ishift+ishift1
2450               ires_old=ires
2451               ibeg=0
2452             else
2453               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2454               ires=ires-ishift+ishift1
2455               ires_old=ires
2456             endif
2457             if (res.eq.'ACE' .or. res.eq.'NHE') then
2458               itype(ires)=10
2459             else
2460               itype(ires)=rescode(ires,res,0)
2461             endif
2462           else
2463             ires=ires-ishift+ishift1
2464           endif
2465 !          write (iout,*) "ires_old",ires_old," ires",ires
2466           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2467 !            ishift1=ishift1+1
2468           endif
2469 !          write (2,*) "ires",ires," res ",res!," ity"!,ity
2470           if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2471              res.eq.'NHE'.and.atom(:2).eq.'HN') then
2472             read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2473 !            write (iout,*) "backbone ",atom
2474 #ifdef DEBUG
2475             write (iout,'(2i3,2x,a,3f8.3)') &
2476             ires,itype(ires),res,(c(j,ires),j=1,3)
2477 #endif
2478             iii=iii+1
2479             do j=1,3
2480               sccor(j,iii)=c(j,ires)
2481             enddo
2482 !            write (*,*) card(23:27),ires,itype(ires)
2483           else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2484                    atom.ne.'N' .and. atom.ne.'C' .and. &
2485                    atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2486                    atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2487 !            write (iout,*) "sidechain ",atom
2488             iii=iii+1
2489             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2490           endif
2491         endif
2492       enddo
2493    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2494       if (ires.eq.0) return
2495 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2496 ! system
2497       nres=ires
2498       do i=2,nres-1
2499 !        write (iout,*) i,itype(i)
2500         if (itype(i).eq.ntyp1) then
2501 !          write (iout,*) "dummy",i,itype(i)
2502           do j=1,3
2503             c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2504 !            c(j,i)=(c(j,i-1)+c(j,i+1))/2
2505             dc(j,i)=c(j,i)
2506           enddo
2507         endif
2508       enddo
2509 ! Calculate the CM of the last side chain.
2510       if (iii.gt.0)  then
2511       if (unres_pdb) then
2512         do j=1,3
2513           dc(j,ires)=sccor(j,iii)
2514         enddo
2515       else
2516         call sccenter(ires,iii,sccor)
2517       endif
2518       endif
2519 !      nres=ires
2520       nsup=nres
2521       nstart_sup=1
2522       if (itype(nres).ne.10) then
2523         nres=nres+1
2524         itype(nres)=ntyp1
2525         if (unres_pdb) then
2526 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2527           call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2528           if (fail) then
2529             e2(1)=0.0d0
2530             e2(2)=1.0d0
2531             e2(3)=0.0d0
2532           endif
2533           do j=1,3
2534             c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
2535           enddo
2536         else
2537         do j=1,3
2538           dcj=c(j,nres-2)-c(j,nres-3)
2539           c(j,nres)=c(j,nres-1)+dcj
2540           c(j,2*nres)=c(j,nres)
2541         enddo
2542         endif
2543       endif
2544 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2545 #ifdef WHAM_RUN
2546       if (nres.ne.nres0) then
2547         write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2548                        " NRES0=",nres0
2549         stop "Error nres value in WHAM input"
2550       endif
2551 #endif
2552 !---------------------------------
2553 !el reallocate tables
2554 !      do i=1,maxres/3
2555 !       do j=1,2
2556 !         hfrag_alloc(j,i)=hfrag(j,i)
2557 !        enddo
2558 !       do j=1,4
2559 !         bfrag_alloc(j,i)=bfrag(j,i)
2560 !        enddo
2561 !      enddo
2562
2563 !      deallocate(hfrag)
2564 !      deallocate(bfrag)
2565 !      allocate(hfrag(2,nres/3)) !(2,maxres/3)
2566 !el      allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2567 !el      allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2568 !      allocate(bfrag(4,nres/3)) !(4,maxres/3)
2569
2570 !      do i=1,nhfrag
2571 !       do j=1,2
2572 !         hfrag(j,i)=hfrag_alloc(j,i)
2573 !        enddo
2574 !      enddo
2575 !      do i=1,nbfrag
2576 !       do j=1,4
2577 !         bfrag(j,i)=bfrag_alloc(j,i)
2578 !        enddo
2579 !      enddo
2580 !el end reallocate tables
2581 !---------------------------------
2582       do i=2,nres-1
2583         do j=1,3
2584           c(j,i+nres)=dc(j,i)
2585         enddo
2586       enddo
2587       do j=1,3
2588         c(j,nres+1)=c(j,1)
2589         c(j,2*nres)=c(j,nres)
2590       enddo
2591       if (itype(1).eq.ntyp1) then
2592         nsup=nsup-1
2593         nstart_sup=2
2594         if (unres_pdb) then
2595 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2596           call refsys(2,3,4,e1,e2,e3,fail)
2597           if (fail) then
2598             e2(1)=0.0d0
2599             e2(2)=1.0d0
2600             e2(3)=0.0d0
2601           endif
2602           do j=1,3
2603             c(j,1)=c(j,2)-3.8d0*e2(j)
2604           enddo
2605         else
2606         do j=1,3
2607           dcj=c(j,4)-c(j,3)
2608           c(j,1)=c(j,2)-dcj
2609           c(j,nres+1)=c(j,1)
2610         enddo
2611         endif
2612       endif
2613 ! Copy the coordinates to reference coordinates
2614 !      do i=1,2*nres
2615 !        do j=1,3
2616 !          cref(j,i)=c(j,i)
2617 !        enddo
2618 !      enddo
2619 ! Calculate internal coordinates.
2620       if (lprn) then
2621       write (iout,'(/a)') &
2622         "Cartesian coordinates of the reference structure"
2623       write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2624        "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2625       do ires=1,nres
2626         write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2627           restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2628           (c(j,ires+nres),j=1,3)
2629       enddo
2630       endif
2631 ! znamy już nres wiec mozna alokowac tablice
2632 ! Calculate internal coordinates.
2633       if(me.eq.king.or..not.out1file)then
2634        write (iout,'(a)') &
2635          "Backbone and SC coordinates as read from the PDB"
2636        do ires=1,nres
2637         write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2638           ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2639           (c(j,nres+ires),j=1,3)
2640        enddo
2641       endif
2642
2643       if(.not.allocated(vbld)) then
2644        allocate(vbld(2*nres))
2645        do i=1,2*nres
2646          vbld(i)=0.d0
2647        enddo
2648       endif
2649       if(.not.allocated(vbld_inv)) then
2650        allocate(vbld_inv(2*nres))
2651        do i=1,2*nres
2652          vbld_inv(i)=0.d0
2653        enddo
2654       endif
2655 !!!el
2656       if(.not.allocated(theta)) then
2657         allocate(theta(nres+2))
2658 !        allocate(phi(nres+2))
2659 !        allocate(alph(nres+2))
2660 !        allocate(omeg(nres+2))
2661         do i=1,nres+2
2662           theta(i)=0.0d0
2663 !          phi(i)=0.0d0
2664 !          alph(i)=0.0d0
2665 !          omeg(i)=0.0d0
2666         enddo
2667       endif
2668 !       allocate(costtab(nres))
2669 !        allocate(sinttab(nres))
2670 !        allocate(cost2tab(nres))
2671 !        allocate(sint2tab(nres))
2672 !        allocate(xxref(nres))
2673 !        allocate(yyref(nres))
2674 !        allocate(zzref(nres)) !(maxres)
2675 !        do i=1,nres
2676 !          costtab(i)=0.0d0
2677 !          sinttab(i)=0.0d0
2678 !          cost2tab(i)=0.0d0
2679 !          sint2tab(i)=0.0d0
2680 !          xxref(i)=0.0d0
2681 !          yyref(i)=0.0d0
2682 !          zzref(i)=0.0d0
2683 !        enddo
2684
2685 !      endif
2686       if(.not.allocated(phi)) allocate(phi(nres+2))
2687       if(.not.allocated(alph)) allocate(alph(nres+2))
2688       if(.not.allocated(omeg)) allocate(omeg(nres+2))
2689       if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2690       if(.not.allocated(phiref)) allocate(phiref(nres+2))
2691       if(.not.allocated(costtab)) allocate(costtab(nres))
2692       if(.not.allocated(sinttab)) allocate(sinttab(nres))
2693       if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2694       if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2695       if(.not.allocated(xxref)) allocate(xxref(nres))
2696       if(.not.allocated(yyref)) allocate(yyref(nres))
2697       if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2698       if(.not.allocated(dc_norm)) then
2699 !      if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2700         allocate(dc_norm(3,0:2*nres+2))
2701         do i=0,2*nres+2
2702           dc_norm(1,i)=0.d0
2703           dc_norm(2,i)=0.d0
2704           dc_norm(3,i)=0.d0
2705         enddo
2706       endif
2707
2708       call int_from_cart(.true.,.false.)
2709       call sc_loc_geom(.true.)
2710 !      call sc_loc_geom(.false.)
2711 ! wczesbiej bylo false
2712       do i=1,nres
2713         thetaref(i)=theta(i)
2714         phiref(i)=phi(i)
2715       enddo
2716 !      do i=1,2*nres
2717 !        vbld_inv(i)=0.d0
2718 !        vbld(i)=0.d0
2719 !      enddo
2720
2721       do i=1,nres-1
2722         do j=1,3
2723           dc(j,i)=c(j,i+1)-c(j,i)
2724           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2725         enddo
2726       enddo
2727       do i=2,nres-1
2728         do j=1,3
2729           dc(j,i+nres)=c(j,i+nres)-c(j,i)
2730           dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2731         enddo
2732 !      write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2733 !        vbld_inv(i+nres)
2734       enddo
2735 !      call chainbuild
2736 ! Copy the coordinates to reference coordinates
2737 ! Splits to single chain if occurs
2738
2739 !      do i=1,2*nres
2740 !        do j=1,3
2741 !          cref(j,i,cou)=c(j,i)
2742 !        enddo
2743 !      enddo
2744 !
2745       if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2746       if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2747       if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2748 !-----------------------------
2749       kkk=1
2750       lll=0
2751       cou=1
2752       do i=1,nres
2753       lll=lll+1
2754 !c      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2755       if (i.gt.1) then
2756       if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2757       chain_length=lll-1
2758       kkk=kkk+1
2759 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2760       lll=1
2761       endif
2762       endif
2763         do j=1,3
2764           cref(j,i,cou)=c(j,i)
2765           cref(j,i+nres,cou)=c(j,i+nres)
2766           if (i.le.nres) then
2767           chain_rep(j,lll,kkk)=c(j,i)
2768           chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2769           endif
2770          enddo
2771       enddo
2772       write (iout,*) chain_length
2773       if (chain_length.eq.0) chain_length=nres
2774       do j=1,3
2775       chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2776       chain_rep(j,chain_length+nres,symetr) &
2777       =chain_rep(j,chain_length+nres,1)
2778       enddo
2779 ! diagnostic
2780 !       write (iout,*) "spraw lancuchy",chain_length,symetr
2781 !       do i=1,4
2782 !         do kkk=1,chain_length
2783 !           write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2784 !         enddo
2785 !        enddo
2786 ! enddiagnostic
2787 ! makes copy of chains
2788         write (iout,*) "symetr", symetr
2789
2790       if (symetr.gt.1) then
2791        call permut(symetr)
2792        nperm=1
2793        do i=1,symetr
2794        nperm=nperm*i
2795        enddo
2796        do i=1,nperm
2797        write(iout,*) (tabperm(i,kkk),kkk=1,4)
2798        enddo
2799        do i=1,nperm
2800         cou=0
2801         do kkk=1,symetr
2802          icha=tabperm(i,kkk)
2803 !         write (iout,*) i,icha
2804          do lll=1,chain_length
2805           cou=cou+1
2806            if (cou.le.nres) then
2807            do j=1,3
2808             kupa=mod(lll,chain_length)
2809             iprzes=(kkk-1)*chain_length+lll
2810             if (kupa.eq.0) kupa=chain_length
2811 !            write (iout,*) "kupa", kupa
2812             cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2813             cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2814           enddo
2815           endif
2816          enddo
2817         enddo
2818        enddo
2819        endif
2820 !-koniec robienia kopii
2821 ! diag
2822       do kkk=1,nperm
2823       write (iout,*) "nowa struktura", nperm
2824       do i=1,nres
2825         write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2826       cref(2,i,kkk),&
2827       cref(3,i,kkk),cref(1,nres+i,kkk),&
2828       cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2829       enddo
2830   100 format (//'              alpha-carbon coordinates       ',&
2831                 '     centroid coordinates'/ &
2832                 '       ', 6X,'X',11X,'Y',11X,'Z', &
2833                                 10X,'X',11X,'Y',11X,'Z')
2834   110 format (a,'(',i3,')',6f12.5)
2835
2836       enddo
2837 !c enddiag
2838       do j=1,nbfrag
2839         do i=1,4
2840          bfrag(i,j)=bfrag(i,j)-ishift
2841         enddo
2842       enddo
2843
2844       do j=1,nhfrag
2845         do i=1,2
2846          hfrag(i,j)=hfrag(i,j)-ishift
2847         enddo
2848       enddo
2849       ishift_pdb=ishift
2850
2851       return
2852       end subroutine readpdb
2853 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2854 !-----------------------------------------------------------------------------
2855 ! readrtns_CSA.F
2856 !-----------------------------------------------------------------------------
2857       subroutine read_control
2858 !
2859 ! Read contorl data
2860 !
2861 !      use geometry_data
2862       use comm_machsw
2863       use energy_data
2864       use control_data
2865       use compare_data
2866       use MCM_data
2867       use map_data
2868       use csa_data
2869       use MD_data
2870       use MPI_data
2871       use random, only: random_init
2872 !      implicit real*8 (a-h,o-z)
2873 !      include 'DIMENSIONS'
2874 #ifdef MP
2875       use prng, only:prng_restart
2876       include 'mpif.h'
2877       logical :: OKRandom!, prng_restart
2878       real(kind=8) :: r1
2879 #endif
2880 !      include 'COMMON.IOUNITS'
2881 !      include 'COMMON.TIME1'
2882 !      include 'COMMON.THREAD'
2883 !      include 'COMMON.SBRIDGE'
2884 !      include 'COMMON.CONTROL'
2885 !      include 'COMMON.MCM'
2886 !      include 'COMMON.MAP'
2887 !      include 'COMMON.HEADER'
2888 !      include 'COMMON.CSA'
2889 !      include 'COMMON.CHAIN'
2890 !      include 'COMMON.MUCA'
2891 !      include 'COMMON.MD'
2892 !      include 'COMMON.FFIELD'
2893 !      include 'COMMON.INTERACT'
2894 !      include 'COMMON.SETUP'
2895 !el      integer :: KDIAG,ICORFL,IXDR
2896 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2897       character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2898         'EVVRSP  ','Givens  ','Jacobi  '/),shape(diagmeth))
2899 !      character(len=80) :: ucase
2900       character(len=640) :: controlcard
2901
2902       real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2903
2904
2905       nglob_csa=0
2906       eglob_csa=1d99
2907       nmin_csa=0
2908       read (INP,'(a)') titel
2909       call card_concat(controlcard,.true.)
2910 !      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2911 !      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2912       call reada(controlcard,'SEED',seed,0.0D0)
2913       call random_init(seed)
2914 ! Set up the time limit (caution! The time must be input in minutes!)
2915       read_cart=index(controlcard,'READ_CART').gt.0
2916       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2917       call readi(controlcard,'SYM',symetr,1)
2918       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2919       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2920       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2921       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2922       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2923       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2924       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2925       call reada(controlcard,'DRMS',drms,0.1D0)
2926       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2927        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2928        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2929        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2930        write (iout,'(a,f10.1)')'DRMS    = ',drms
2931        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2932        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2933       endif
2934       call readi(controlcard,'NZ_START',nz_start,0)
2935       call readi(controlcard,'NZ_END',nz_end,0)
2936       call readi(controlcard,'IZ_SC',iz_sc,0)
2937       timlim=60.0D0*timlim
2938       safety = 60.0d0*safety
2939       timem=timlim
2940       modecalc=0
2941       call reada(controlcard,"T_BATH",t_bath,300.0d0)
2942       minim=(index(controlcard,'MINIMIZE').gt.0)
2943       dccart=(index(controlcard,'CART').gt.0)
2944       overlapsc=(index(controlcard,'OVERLAP').gt.0)
2945       overlapsc=.not.overlapsc
2946       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
2947       searchsc=.not.searchsc
2948       sideadd=(index(controlcard,'SIDEADD').gt.0)
2949       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
2950       outpdb=(index(controlcard,'PDBOUT').gt.0)
2951       outmol2=(index(controlcard,'MOL2OUT').gt.0)
2952       pdbref=(index(controlcard,'PDBREF').gt.0)
2953       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
2954       indpdb=index(controlcard,'PDBSTART')
2955       extconf=(index(controlcard,'EXTCONF').gt.0)
2956       call readi(controlcard,'IPRINT',iprint,0)
2957       call readi(controlcard,'MAXGEN',maxgen,10000)
2958       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
2959       call readi(controlcard,"KDIAG",kdiag,0)
2960       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
2961       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
2962        write (iout,*) "RESCALE_MODE",rescale_mode
2963       split_ene=index(controlcard,'SPLIT_ENE').gt.0
2964       if (index(controlcard,'REGULAR').gt.0.0D0) then
2965         call reada(controlcard,'WEIDIS',weidis,0.1D0)
2966         modecalc=1
2967         refstr=.true.
2968       endif
2969       if (index(controlcard,'CHECKGRAD').gt.0) then
2970         modecalc=5
2971         if (index(controlcard,'CART').gt.0) then
2972           icheckgrad=1
2973         elseif (index(controlcard,'CARINT').gt.0) then
2974           icheckgrad=2
2975         else
2976           icheckgrad=3
2977         endif
2978       elseif (index(controlcard,'THREAD').gt.0) then
2979         modecalc=2
2980         call readi(controlcard,'THREAD',nthread,0)
2981         if (nthread.gt.0) then
2982           call reada(controlcard,'WEIDIS',weidis,0.1D0)
2983         else
2984           if (fg_rank.eq.0) &
2985           write (iout,'(a)')'A number has to follow the THREAD keyword.'
2986           stop 'Error termination in Read_Control.'
2987         endif
2988       else if (index(controlcard,'MCMA').gt.0) then
2989         modecalc=3
2990       else if (index(controlcard,'MCEE').gt.0) then
2991         modecalc=6
2992       else if (index(controlcard,'MULTCONF').gt.0) then
2993         modecalc=4
2994       else if (index(controlcard,'MAP').gt.0) then
2995         modecalc=7
2996         call readi(controlcard,'MAP',nmap,0)
2997       else if (index(controlcard,'CSA').gt.0) then
2998         modecalc=8
2999 !rc      else if (index(controlcard,'ZSCORE').gt.0) then
3000 !rc
3001 !rc  ZSCORE is rm from UNRES, modecalc=9 is available
3002 !rc
3003 !rc        modecalc=9
3004 !fcm      else if (index(controlcard,'MCMF').gt.0) then
3005 !fmc        modecalc=10
3006       else if (index(controlcard,'SOFTREG').gt.0) then
3007         modecalc=11
3008       else if (index(controlcard,'CHECK_BOND').gt.0) then
3009         modecalc=-1
3010       else if (index(controlcard,'TEST').gt.0) then
3011         modecalc=-2
3012       else if (index(controlcard,'MD').gt.0) then
3013         modecalc=12
3014       else if (index(controlcard,'RE ').gt.0) then
3015         modecalc=14
3016       endif
3017
3018       lmuca=index(controlcard,'MUCA').gt.0
3019       call readi(controlcard,'MUCADYN',mucadyn,0)
3020       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3021       if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3022        then
3023        write (iout,*) 'MUCADYN=',mucadyn
3024        write (iout,*) 'MUCASMOOTH=',muca_smooth
3025       endif
3026
3027       iscode=index(controlcard,'ONE_LETTER')
3028       indphi=index(controlcard,'PHI')
3029       indback=index(controlcard,'BACK')
3030       iranconf=index(controlcard,'RAND_CONF')
3031       i2ndstr=index(controlcard,'USE_SEC_PRED')
3032       gradout=index(controlcard,'GRADOUT').gt.0
3033       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3034       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3035       if (me.eq.king .or. .not.out1file ) &
3036         write (iout,*) "DISTCHAINMAX",distchainmax
3037
3038       if(me.eq.king.or..not.out1file) &
3039        write (iout,'(2a)') diagmeth(kdiag),&
3040         ' routine used to diagonalize matrices.'
3041       return
3042       end subroutine read_control
3043 !-----------------------------------------------------------------------------
3044       subroutine read_REMDpar
3045 !
3046 ! Read REMD settings
3047 !
3048 !       use control
3049 !       use energy
3050 !       use geometry
3051       use REMD_data
3052       use MPI_data
3053       use control_data, only:out1file
3054 !      implicit real*8 (a-h,o-z)
3055 !      include 'DIMENSIONS'
3056 !      include 'COMMON.IOUNITS'
3057 !      include 'COMMON.TIME1'
3058 !      include 'COMMON.MD'
3059       use MD_data
3060 !el #ifndef LANG0
3061 !el      include 'COMMON.LANGEVIN'
3062 !el #else
3063 !el      include 'COMMON.LANGEVIN.lang0'
3064 !el #endif
3065 !      include 'COMMON.INTERACT'
3066 !      include 'COMMON.NAMES'
3067 !      include 'COMMON.GEO'
3068 !      include 'COMMON.REMD'
3069 !      include 'COMMON.CONTROL'
3070 !      include 'COMMON.SETUP'
3071 !      character(len=80) :: ucase
3072       character(len=320) :: controlcard
3073       character(len=3200) :: controlcard1
3074       integer :: iremd_m_total
3075 !el local variables
3076       integer :: i
3077 !     real(kind=8) :: var,ene
3078
3079       if(me.eq.king.or..not.out1file) &
3080        write (iout,*) "REMD setup"
3081
3082       call card_concat(controlcard,.true.)
3083       call readi(controlcard,"NREP",nrep,3)
3084       call readi(controlcard,"NSTEX",nstex,1000)
3085       call reada(controlcard,"RETMIN",retmin,10.0d0)
3086       call reada(controlcard,"RETMAX",retmax,1000.0d0)
3087       mremdsync=(index(controlcard,'SYNC').gt.0)
3088       call readi(controlcard,"NSYN",i_sync_step,100)
3089       restart1file=(index(controlcard,'REST1FILE').gt.0)
3090       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3091       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3092       if(max_cache_traj_use.gt.max_cache_traj) &
3093                 max_cache_traj_use=max_cache_traj
3094       if(me.eq.king.or..not.out1file) then
3095 !d       if (traj1file) then
3096 !rc caching is in testing - NTWX is not ignored
3097 !d        write (iout,*) "NTWX value is ignored"
3098 !d        write (iout,*) "  trajectory is stored to one file by master"
3099 !d        write (iout,*) "  before exchange at NSTEX intervals"
3100 !d       endif
3101        write (iout,*) "NREP= ",nrep
3102        write (iout,*) "NSTEX= ",nstex
3103        write (iout,*) "SYNC= ",mremdsync
3104        write (iout,*) "NSYN= ",i_sync_step
3105        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3106       endif
3107       remd_tlist=.false.
3108       allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3109       if (index(controlcard,'TLIST').gt.0) then
3110          remd_tlist=.true.
3111          call card_concat(controlcard1,.true.)
3112          read(controlcard1,*) (remd_t(i),i=1,nrep)
3113          if(me.eq.king.or..not.out1file) &
3114           write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3115       endif
3116       remd_mlist=.false.
3117       if (index(controlcard,'MLIST').gt.0) then
3118          remd_mlist=.true.
3119          call card_concat(controlcard1,.true.)
3120          read(controlcard1,*) (remd_m(i),i=1,nrep)
3121          if(me.eq.king.or..not.out1file) then
3122           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3123           iremd_m_total=0
3124           do i=1,nrep
3125            iremd_m_total=iremd_m_total+remd_m(i)
3126           enddo
3127           write (iout,*) 'Total number of replicas ',iremd_m_total
3128          endif
3129       endif
3130       if(me.eq.king.or..not.out1file) &
3131        write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3132       return
3133       end subroutine read_REMDpar
3134 !-----------------------------------------------------------------------------
3135       subroutine read_MDpar
3136 !
3137 ! Read MD settings
3138 !
3139       use control_data, only: r_cut,rlamb,out1file
3140       use energy_data
3141       use geometry_data, only: pi
3142       use MPI_data
3143 !      implicit real*8 (a-h,o-z)
3144 !      include 'DIMENSIONS'
3145 !      include 'COMMON.IOUNITS'
3146 !      include 'COMMON.TIME1'
3147 !      include 'COMMON.MD'
3148       use MD_data
3149 !el #ifndef LANG0
3150 !el      include 'COMMON.LANGEVIN'
3151 !el #else
3152 !el      include 'COMMON.LANGEVIN.lang0'
3153 !el #endif
3154 !      include 'COMMON.INTERACT'
3155 !      include 'COMMON.NAMES'
3156 !      include 'COMMON.GEO'
3157 !      include 'COMMON.SETUP'
3158 !      include 'COMMON.CONTROL'
3159 !      include 'COMMON.SPLITELE'
3160 !      character(len=80) :: ucase
3161       character(len=320) :: controlcard
3162 !el local variables
3163       integer :: i
3164       real(kind=8) :: eta
3165
3166       call card_concat(controlcard,.true.)
3167       call readi(controlcard,"NSTEP",n_timestep,1000000)
3168       call readi(controlcard,"NTWE",ntwe,100)
3169       call readi(controlcard,"NTWX",ntwx,1000)
3170       call reada(controlcard,"DT",d_time,1.0d-1)
3171       call reada(controlcard,"DVMAX",dvmax,2.0d1)
3172       call reada(controlcard,"DAMAX",damax,1.0d1)
3173       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3174       call readi(controlcard,"LANG",lang,0)
3175       RESPA = index(controlcard,"RESPA") .gt. 0
3176       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3177       ntime_split0=ntime_split
3178       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3179       ntime_split0=ntime_split
3180       call reada(controlcard,"R_CUT",r_cut,2.0d0)
3181       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3182       rest = index(controlcard,"REST").gt.0
3183       tbf = index(controlcard,"TBF").gt.0
3184       usampl = index(controlcard,"USAMPL").gt.0
3185       mdpdb = index(controlcard,"MDPDB").gt.0
3186       call reada(controlcard,"T_BATH",t_bath,300.0d0)
3187       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3188       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3189       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3190       if (count_reset_moment.eq.0) count_reset_moment=1000000000
3191       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3192       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3193       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3194       if (count_reset_vel.eq.0) count_reset_vel=1000000000
3195       large = index(controlcard,"LARGE").gt.0
3196       print_compon = index(controlcard,"PRINT_COMPON").gt.0
3197       rattle = index(controlcard,"RATTLE").gt.0
3198 !  if performing umbrella sampling, fragments constrained are read from the fragment file
3199       nset=0
3200       if(usampl) then
3201         call read_fragments
3202       endif
3203
3204       if(me.eq.king.or..not.out1file) then
3205        write (iout,*)
3206        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3207        write (iout,*)
3208        write (iout,'(a)') "The units are:"
3209        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3210        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3211         " acceleration: angstrom/(48.9 fs)**2"
3212        write (iout,'(a)') "energy: kcal/mol, temperature: K"
3213        write (iout,*)
3214        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3215        write (iout,'(a60,f10.5,a)') &
3216         "Initial time step of numerical integration:",d_time,&
3217         " natural units"
3218        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3219        if (RESPA) then
3220         write (iout,'(2a,i4,a)') &
3221           "A-MTS algorithm used; initial time step for fast-varying",&
3222           " short-range forces split into",ntime_split," steps."
3223         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3224          r_cut," lambda",rlamb
3225        endif
3226        write (iout,'(2a,f10.5)') &
3227         "Maximum acceleration threshold to reduce the time step",&
3228         "/increase split number:",damax
3229        write (iout,'(2a,f10.5)') &
3230         "Maximum predicted energy drift to reduce the timestep",&
3231         "/increase split number:",edriftmax
3232        write (iout,'(a60,f10.5)') &
3233        "Maximum velocity threshold to reduce velocities:",dvmax
3234        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3235        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3236        if (rattle) write (iout,'(a60)') &
3237         "Rattle algorithm used to constrain the virtual bonds"
3238       endif
3239       reset_fricmat=1000
3240       if (lang.gt.0) then
3241         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3242         call reada(controlcard,"RWAT",rwat,1.4d0)
3243         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3244         surfarea=index(controlcard,"SURFAREA").gt.0
3245         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3246         if(me.eq.king.or..not.out1file)then
3247          write (iout,'(/a,$)') "Langevin dynamics calculation"
3248          if (lang.eq.1) then
3249           write (iout,'(a/)') &
3250             " with direct integration of Langevin equations"
3251          else if (lang.eq.2) then
3252           write (iout,'(a/)') " with TINKER stochasic MD integrator"
3253          else if (lang.eq.3) then
3254           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3255          else if (lang.eq.4) then
3256           write (iout,'(a/)') " in overdamped mode"
3257          else
3258           write (iout,'(//a,i5)') &
3259             "=========== ERROR: Unknown Langevin dynamics mode:",lang
3260           stop
3261          endif
3262          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3263          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3264          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3265          write (iout,'(a60,f10.5)') &
3266          "Scaling factor of the friction forces:",scal_fric
3267          if (surfarea) write (iout,'(2a,i10,a)') &
3268            "Friction coefficients will be scaled by solvent-accessible",&
3269            " surface area every",reset_fricmat," steps."
3270         endif
3271 ! Calculate friction coefficients and bounds of stochastic forces
3272         eta=6*pi*cPoise*etawat
3273         if(me.eq.king.or..not.out1file) &
3274          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3275           eta
3276         gamp=scal_fric*(pstok+rwat)*eta
3277         stdfp=dsqrt(2*Rb*t_bath/d_time)
3278         allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3279         do i=1,ntyp
3280           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3281           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3282         enddo
3283         if(me.eq.king.or..not.out1file)then
3284          write (iout,'(/2a/)') &
3285          "Radii of site types and friction coefficients and std's of",&
3286          " stochastic forces of fully exposed sites"
3287          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3288          do i=1,ntyp
3289           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3290            gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3291          enddo
3292         endif
3293       else if (tbf) then
3294         if(me.eq.king.or..not.out1file)then
3295          write (iout,'(a)') "Berendsen bath calculation"
3296          write (iout,'(a60,f10.5)') "Temperature:",t_bath
3297          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3298          if (reset_moment) &
3299          write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3300          count_reset_moment," steps"
3301          if (reset_vel) &
3302           write (iout,'(a,i10,a)') &
3303           "Velocities will be reset at random every",count_reset_vel,&
3304          " steps"
3305         endif
3306       else
3307         if(me.eq.king.or..not.out1file) &
3308          write (iout,'(a31)') "Microcanonical mode calculation"
3309       endif
3310       if(me.eq.king.or..not.out1file)then
3311        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3312        if (usampl) then
3313           write(iout,*) "MD running with constraints."
3314           write(iout,*) "Equilibration time ", eq_time, " mtus."
3315           write(iout,*) "Constraining ", nfrag," fragments."
3316           write(iout,*) "Length of each fragment, weight and q0:"
3317           do iset=1,nset
3318            write (iout,*) "Set of restraints #",iset
3319            do i=1,nfrag
3320               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3321                  ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3322            enddo
3323            write(iout,*) "constraints between ", npair, "fragments."
3324            write(iout,*) "constraint pairs, weights and q0:"
3325            do i=1,npair
3326             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3327                    ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3328            enddo
3329            write(iout,*) "angle constraints within ", nfrag_back,&
3330             "backbone fragments."
3331            write(iout,*) "fragment, weights:"
3332            do i=1,nfrag_back
3333             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3334                ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3335                wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3336            enddo
3337           enddo
3338         iset=mod(kolor,nset)+1
3339        endif
3340       endif
3341       if(me.eq.king.or..not.out1file) &
3342        write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3343       return
3344       end subroutine read_MDpar
3345 !-----------------------------------------------------------------------------
3346       subroutine map_read
3347
3348       use map_data
3349 !      implicit real*8 (a-h,o-z)
3350 !      include 'DIMENSIONS'
3351 !      include 'COMMON.MAP'
3352 !      include 'COMMON.IOUNITS'
3353       character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3354       character(len=80) :: mapcard      !,ucase
3355 !el local variables
3356       integer :: imap
3357 !     real(kind=8) :: var,ene
3358
3359       do imap=1,nmap
3360         read (inp,'(a)') mapcard
3361         mapcard=ucase(mapcard)
3362         if (index(mapcard,'PHI').gt.0) then
3363           kang(imap)=1
3364         else if (index(mapcard,'THE').gt.0) then
3365           kang(imap)=2
3366         else if (index(mapcard,'ALP').gt.0) then
3367           kang(imap)=3
3368         else if (index(mapcard,'OME').gt.0) then
3369           kang(imap)=4
3370         else
3371           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3372           stop 'Error - illegal variable spec in MAP card.'
3373         endif
3374         call readi (mapcard,'RES1',res1(imap),0)
3375         call readi (mapcard,'RES2',res2(imap),0)
3376         if (res1(imap).eq.0) then
3377           res1(imap)=res2(imap)
3378         else if (res2(imap).eq.0) then
3379           res2(imap)=res1(imap)
3380         endif
3381         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3382           write (iout,'(a)') &
3383           'Error - illegal definition of variable group in MAP.'
3384           stop 'Error - illegal definition of variable group in MAP.'
3385         endif
3386         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3387         call reada(mapcard,'TO',ang_to(imap),0.0D0)
3388         call readi(mapcard,'NSTEP',nstep(imap),0)
3389         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3390           write (iout,'(a)') &
3391            'Illegal boundary and/or step size specification in MAP.'
3392           stop 'Illegal boundary and/or step size specification in MAP.'
3393         endif
3394       enddo ! imap
3395       return
3396       end subroutine map_read
3397 !-----------------------------------------------------------------------------
3398       subroutine csaread
3399
3400       use control_data, only: vdisulf
3401       use csa_data
3402 !      implicit real*8 (a-h,o-z)
3403 !      include 'DIMENSIONS'
3404 !      include 'COMMON.IOUNITS'
3405 !      include 'COMMON.GEO'
3406 !      include 'COMMON.CSA'
3407 !      include 'COMMON.BANK'
3408 !      include 'COMMON.CONTROL'
3409 !      character(len=80) :: ucase
3410       character(len=620) :: mcmcard
3411 !el local variables
3412 !     integer :: ntf,ik,iw_pdb
3413 !     real(kind=8) :: var,ene
3414
3415       call card_concat(mcmcard,.true.)
3416
3417       call readi(mcmcard,'NCONF',nconf,50)
3418       call readi(mcmcard,'NADD',nadd,0)
3419       call readi(mcmcard,'JSTART',jstart,1)
3420       call readi(mcmcard,'JEND',jend,1)
3421       call readi(mcmcard,'NSTMAX',nstmax,500000)
3422       call readi(mcmcard,'N0',n0,1)
3423       call readi(mcmcard,'N1',n1,6)
3424       call readi(mcmcard,'N2',n2,4)
3425       call readi(mcmcard,'N3',n3,0)
3426       call readi(mcmcard,'N4',n4,0)
3427       call readi(mcmcard,'N5',n5,0)
3428       call readi(mcmcard,'N6',n6,10)
3429       call readi(mcmcard,'N7',n7,0)
3430       call readi(mcmcard,'N8',n8,0)
3431       call readi(mcmcard,'N9',n9,0)
3432       call readi(mcmcard,'N14',n14,0)
3433       call readi(mcmcard,'N15',n15,0)
3434       call readi(mcmcard,'N16',n16,0)
3435       call readi(mcmcard,'N17',n17,0)
3436       call readi(mcmcard,'N18',n18,0)
3437
3438       vdisulf=(index(mcmcard,'DYNSS').gt.0)
3439
3440       call readi(mcmcard,'NDIFF',ndiff,2)
3441       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3442       call readi(mcmcard,'IS1',is1,1)
3443       call readi(mcmcard,'IS2',is2,8)
3444       call readi(mcmcard,'NRAN0',nran0,4)
3445       call readi(mcmcard,'NRAN1',nran1,2)
3446       call readi(mcmcard,'IRR',irr,1)
3447       call readi(mcmcard,'NSEED',nseed,20)
3448       call readi(mcmcard,'NTOTAL',ntotal,10000)
3449       call reada(mcmcard,'CUT1',cut1,2.0d0)
3450       call reada(mcmcard,'CUT2',cut2,5.0d0)
3451       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3452       call readi(mcmcard,'ICMAX',icmax,3)
3453       call readi(mcmcard,'IRESTART',irestart,0)
3454 !!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
3455       ntbankm=0
3456 !!bankt
3457       call reada(mcmcard,'DELE',dele,20.0d0)
3458       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3459       call readi(mcmcard,'IREF',iref,0)
3460       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3461       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3462       call readi(mcmcard,'NCONF_IN',nconf_in,0)
3463       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3464       write (iout,*) "NCONF_IN",nconf_in
3465       return
3466       end subroutine csaread
3467 !-----------------------------------------------------------------------------
3468       subroutine mcmread
3469
3470       use mcm_data
3471       use control_data, only: MaxMoveType
3472       use MD_data
3473       use minim_data
3474 !      implicit real*8 (a-h,o-z)
3475 !      include 'DIMENSIONS'
3476 !      include 'COMMON.MCM'
3477 !      include 'COMMON.MCE'
3478 !      include 'COMMON.IOUNITS'
3479 !      character(len=80) :: ucase
3480       character(len=320) :: mcmcard
3481 !el local variables
3482       integer :: i
3483 !     real(kind=8) :: var,ene
3484
3485       call card_concat(mcmcard,.true.)
3486       call readi(mcmcard,'MAXACC',maxacc,100)
3487       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3488       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3489       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3490       call readi(mcmcard,'MAXREPM',maxrepm,200)
3491       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3492       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3493       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3494       call reada(mcmcard,'E_UP',e_up,5.0D0)
3495       call reada(mcmcard,'DELTE',delte,0.1D0)
3496       call readi(mcmcard,'NSWEEP',nsweep,5)
3497       call readi(mcmcard,'NSTEPH',nsteph,0)
3498       call readi(mcmcard,'NSTEPC',nstepc,0)
3499       call reada(mcmcard,'TMIN',tmin,298.0D0)
3500       call reada(mcmcard,'TMAX',tmax,298.0D0)
3501       call readi(mcmcard,'NWINDOW',nwindow,0)
3502       call readi(mcmcard,'PRINT_MC',print_mc,0)
3503       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3504       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3505       ent_read=(index(mcmcard,'ENT_READ').gt.0)
3506       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3507       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3508       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3509       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3510       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3511       if (nwindow.gt.0) then
3512         allocate(winstart(nwindow))     !!el (maxres)
3513         allocate(winend(nwindow))       !!el
3514         allocate(winlen(nwindow))       !!el
3515         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3516         do i=1,nwindow
3517           winlen(i)=winend(i)-winstart(i)+1
3518         enddo
3519       endif
3520       if (tmax.lt.tmin) tmax=tmin
3521       if (tmax.eq.tmin) then
3522         nstepc=0
3523         nsteph=0
3524       endif
3525       if (nstepc.gt.0 .and. nsteph.gt.0) then
3526         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3527         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3528       endif
3529       allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3530 ! Probabilities of different move types
3531       sumpro_type(0)=0.0D0
3532       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3533       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3534       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3535       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
3536       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3537       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3538       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3539       do i=1,MaxMoveType
3540         print *,'i',i,' sumprotype',sumpro_type(i)
3541         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3542         print *,'i',i,' sumprotype',sumpro_type(i)
3543       enddo
3544       return
3545       end subroutine mcmread
3546 !-----------------------------------------------------------------------------
3547       subroutine read_minim
3548
3549       use minim_data
3550 !      implicit real*8 (a-h,o-z)
3551 !      include 'DIMENSIONS'
3552 !      include 'COMMON.MINIM'
3553 !      include 'COMMON.IOUNITS'
3554 !      character(len=80) :: ucase
3555       character(len=320) :: minimcard
3556 !el local variables
3557 !     integer :: ntf,ik,iw_pdb
3558 !     real(kind=8) :: var,ene
3559
3560       call card_concat(minimcard,.true.)
3561       call readi(minimcard,'MAXMIN',maxmin,2000)
3562       call readi(minimcard,'MAXFUN',maxfun,5000)
3563       call readi(minimcard,'MINMIN',minmin,maxmin)
3564       call readi(minimcard,'MINFUN',minfun,maxmin)
3565       call reada(minimcard,'TOLF',tolf,1.0D-2)
3566       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3567       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3568       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3569       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3570       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3571                'Options in energy minimization:'
3572       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3573        'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3574        'MinMin:',MinMin,' MinFun:',MinFun,&
3575        ' TolF:',TolF,' RTolF:',RTolF
3576       return
3577       end subroutine read_minim
3578 !-----------------------------------------------------------------------------
3579       subroutine openunits
3580
3581       use energy_data, only: usampl
3582       use csa_data
3583       use MPI_data
3584       use control_data, only:out1file
3585       use control, only: getenv_loc
3586 !      implicit real*8 (a-h,o-z)
3587 !      include 'DIMENSIONS'
3588 #ifdef MPI
3589       include 'mpif.h'
3590       character(len=16) :: form,nodename
3591       integer :: nodelen,ierror,npos
3592 #endif
3593 !      include 'COMMON.SETUP'
3594 !      include 'COMMON.IOUNITS'
3595 !      include 'COMMON.MD'
3596 !      include 'COMMON.CONTROL'
3597       integer :: lenpre,lenpot,lentmp   !,ilen
3598 !el      external ilen
3599       character(len=3) :: out1file_text !,ucase
3600       character(len=3) :: ll
3601 !el      external ucase
3602 !el local variables
3603 !     integer :: ntf,ik,iw_pdb
3604 !     real(kind=8) :: var,ene
3605 !
3606 !      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3607       call getenv_loc("PREFIX",prefix)
3608       pref_orig = prefix
3609       call getenv_loc("POT",pot)
3610       call getenv_loc("DIRTMP",tmpdir)
3611       call getenv_loc("CURDIR",curdir)
3612       call getenv_loc("OUT1FILE",out1file_text)
3613 !      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3614       out1file_text=ucase(out1file_text)
3615       if (out1file_text(1:1).eq."Y") then
3616         out1file=.true.
3617       else
3618         out1file=fg_rank.gt.0
3619       endif
3620       lenpre=ilen(prefix)
3621       lenpot=ilen(pot)
3622       lentmp=ilen(tmpdir)
3623       if (lentmp.gt.0) then
3624           write (*,'(80(1h!))')
3625           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
3626           write (*,'(80(1h!))')
3627           write (*,*)"All output files will be on node /tmp directory."
3628 #ifdef MPI
3629         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3630         if (me.eq.king) then
3631           write (*,*) "The master node is ",nodename
3632         else if (fg_rank.eq.0) then
3633           write (*,*) "I am the CG slave node ",nodename
3634         else
3635           write (*,*) "I am the FG slave node ",nodename
3636         endif
3637 #endif
3638         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3639         lenpre = lentmp+lenpre+1
3640       endif
3641       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3642 ! Get the names and open the input files
3643 #if defined(WINIFL) || defined(WINPGI)
3644       open(1,file=pref_orig(:ilen(pref_orig))// &
3645         '.inp',status='old',readonly,shared)
3646        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3647 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3648 ! Get parameter filenames and open the parameter files.
3649       call getenv_loc('BONDPAR',bondname)
3650       open (ibond,file=bondname,status='old',readonly,shared)
3651       call getenv_loc('THETPAR',thetname)
3652       open (ithep,file=thetname,status='old',readonly,shared)
3653       call getenv_loc('ROTPAR',rotname)
3654       open (irotam,file=rotname,status='old',readonly,shared)
3655       call getenv_loc('TORPAR',torname)
3656       open (itorp,file=torname,status='old',readonly,shared)
3657       call getenv_loc('TORDPAR',tordname)
3658       open (itordp,file=tordname,status='old',readonly,shared)
3659       call getenv_loc('FOURIER',fouriername)
3660       open (ifourier,file=fouriername,status='old',readonly,shared)
3661       call getenv_loc('ELEPAR',elename)
3662       open (ielep,file=elename,status='old',readonly,shared)
3663       call getenv_loc('SIDEPAR',sidename)
3664       open (isidep,file=sidename,status='old',readonly,shared)
3665 #elif (defined CRAY) || (defined AIX)
3666       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3667         action='read')
3668 !      print *,"Processor",myrank," opened file 1"
3669       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3670 !      print *,"Processor",myrank," opened file 9"
3671 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3672 ! Get parameter filenames and open the parameter files.
3673       call getenv_loc('BONDPAR',bondname)
3674       open (ibond,file=bondname,status='old',action='read')
3675 !      print *,"Processor",myrank," opened file IBOND"
3676       call getenv_loc('THETPAR',thetname)
3677       open (ithep,file=thetname,status='old',action='read')
3678 !      print *,"Processor",myrank," opened file ITHEP"
3679       call getenv_loc('ROTPAR',rotname)
3680       open (irotam,file=rotname,status='old',action='read')
3681 !      print *,"Processor",myrank," opened file IROTAM"
3682       call getenv_loc('TORPAR',torname)
3683       open (itorp,file=torname,status='old',action='read')
3684 !      print *,"Processor",myrank," opened file ITORP"
3685       call getenv_loc('TORDPAR',tordname)
3686       open (itordp,file=tordname,status='old',action='read')
3687 !      print *,"Processor",myrank," opened file ITORDP"
3688       call getenv_loc('SCCORPAR',sccorname)
3689       open (isccor,file=sccorname,status='old',action='read')
3690 !      print *,"Processor",myrank," opened file ISCCOR"
3691       call getenv_loc('FOURIER',fouriername)
3692       open (ifourier,file=fouriername,status='old',action='read')
3693 !      print *,"Processor",myrank," opened file IFOURIER"
3694       call getenv_loc('ELEPAR',elename)
3695       open (ielep,file=elename,status='old',action='read')
3696 !      print *,"Processor",myrank," opened file IELEP"
3697       call getenv_loc('SIDEPAR',sidename)
3698       open (isidep,file=sidename,status='old',action='read')
3699 !      print *,"Processor",myrank," opened file ISIDEP"
3700 !      print *,"Processor",myrank," opened parameter files"
3701 #elif (defined G77)
3702       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3703       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3704 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3705 ! Get parameter filenames and open the parameter files.
3706       call getenv_loc('BONDPAR',bondname)
3707       open (ibond,file=bondname,status='old')
3708       call getenv_loc('THETPAR',thetname)
3709       open (ithep,file=thetname,status='old')
3710       call getenv_loc('ROTPAR',rotname)
3711       open (irotam,file=rotname,status='old')
3712       call getenv_loc('TORPAR',torname)
3713       open (itorp,file=torname,status='old')
3714       call getenv_loc('TORDPAR',tordname)
3715       open (itordp,file=tordname,status='old')
3716       call getenv_loc('SCCORPAR',sccorname)
3717       open (isccor,file=sccorname,status='old')
3718       call getenv_loc('FOURIER',fouriername)
3719       open (ifourier,file=fouriername,status='old')
3720       call getenv_loc('ELEPAR',elename)
3721       open (ielep,file=elename,status='old')
3722       call getenv_loc('SIDEPAR',sidename)
3723       open (isidep,file=sidename,status='old')
3724 #else
3725       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3726         readonly)
3727        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3728 !      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3729 ! Get parameter filenames and open the parameter files.
3730       call getenv_loc('BONDPAR',bondname)
3731       open (ibond,file=bondname,status='old',action='read')
3732       call getenv_loc('THETPAR',thetname)
3733       open (ithep,file=thetname,status='old',action='read')
3734       call getenv_loc('ROTPAR',rotname)
3735       open (irotam,file=rotname,status='old',action='read')
3736       call getenv_loc('TORPAR',torname)
3737       open (itorp,file=torname,status='old',action='read')
3738       call getenv_loc('TORDPAR',tordname)
3739       open (itordp,file=tordname,status='old',action='read')
3740       call getenv_loc('SCCORPAR',sccorname)
3741       open (isccor,file=sccorname,status='old',action='read')
3742 #ifndef CRYST_THETA
3743       call getenv_loc('THETPARPDB',thetname_pdb)
3744       print *,"thetname_pdb ",thetname_pdb
3745       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3746       print *,ithep_pdb," opened"
3747 #endif
3748       call getenv_loc('FOURIER',fouriername)
3749       open (ifourier,file=fouriername,status='old',readonly)
3750       call getenv_loc('ELEPAR',elename)
3751       open (ielep,file=elename,status='old',readonly)
3752       call getenv_loc('SIDEPAR',sidename)
3753       open (isidep,file=sidename,status='old',readonly)
3754 #ifndef CRYST_SC
3755       call getenv_loc('ROTPARPDB',rotname_pdb)
3756       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3757 #endif
3758 #endif
3759 #ifndef OLDSCP
3760 !
3761 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3762 ! Use -DOLDSCP to use hard-coded constants instead.
3763 !
3764       call getenv_loc('SCPPAR',scpname)
3765 #if defined(WINIFL) || defined(WINPGI)
3766       open (iscpp,file=scpname,status='old',readonly,shared)
3767 #elif (defined CRAY)  || (defined AIX)
3768       open (iscpp,file=scpname,status='old',action='read')
3769 #elif (defined G77)
3770       open (iscpp,file=scpname,status='old')
3771 #else
3772       open (iscpp,file=scpname,status='old',action='read')
3773 #endif
3774 #endif
3775       call getenv_loc('PATTERN',patname)
3776 #if defined(WINIFL) || defined(WINPGI)
3777       open (icbase,file=patname,status='old',readonly,shared)
3778 #elif (defined CRAY)  || (defined AIX)
3779       open (icbase,file=patname,status='old',action='read')
3780 #elif (defined G77)
3781       open (icbase,file=patname,status='old')
3782 #else
3783       open (icbase,file=patname,status='old',action='read')
3784 #endif
3785 #ifdef MPI
3786 ! Open output file only for CG processes
3787 !      print *,"Processor",myrank," fg_rank",fg_rank
3788       if (fg_rank.eq.0) then
3789
3790       if (nodes.eq.1) then
3791         npos=3
3792       else
3793         npos = dlog10(dfloat(nodes-1))+1
3794       endif
3795       if (npos.lt.3) npos=3
3796       write (liczba,'(i1)') npos
3797       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3798         //')'
3799       write (liczba,form) me
3800       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3801         liczba(:ilen(liczba))
3802       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3803         //'.int'
3804       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3805         //'.pdb'
3806       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3807         liczba(:ilen(liczba))//'.mol2'
3808       statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3809         liczba(:ilen(liczba))//'.stat'
3810       if (lentmp.gt.0) &
3811         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3812             //liczba(:ilen(liczba))//'.stat')
3813       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3814         //'.rst'
3815       if(usampl) then
3816           qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3817        liczba(:ilen(liczba))//'.const'
3818       endif
3819
3820       endif
3821 #else
3822       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3823       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3824       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3825       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3826       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3827       if (lentmp.gt.0) &
3828         call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3829           '.stat')
3830       rest2name=prefix(:ilen(prefix))//'.rst'
3831       if(usampl) then
3832          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3833       endif
3834 #endif
3835 #if defined(AIX) || defined(PGI)
3836       if (me.eq.king .or. .not. out1file) &
3837          open(iout,file=outname,status='unknown')
3838 #ifdef DEBUG
3839       if (fg_rank.gt.0) then
3840         write (liczba,'(i3.3)') myrank/nfgtasks
3841         write (ll,'(bz,i3.3)') fg_rank
3842         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3843          status='unknown')
3844       endif
3845 #endif
3846       if(me.eq.king) then
3847        open(igeom,file=intname,status='unknown',position='append')
3848        open(ipdb,file=pdbname,status='unknown')
3849        open(imol2,file=mol2name,status='unknown')
3850        open(istat,file=statname,status='unknown',position='append')
3851       else
3852 !1out       open(iout,file=outname,status='unknown')
3853       endif
3854 #else
3855       if (me.eq.king .or. .not.out1file) &
3856           open(iout,file=outname,status='unknown')
3857 #ifdef DEBUG
3858       if (fg_rank.gt.0) then
3859         write (liczba,'(i3.3)') myrank/nfgtasks
3860         write (ll,'(bz,i3.3)') fg_rank
3861         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3862          status='unknown')
3863       endif
3864 #endif
3865       if(me.eq.king) then
3866        open(igeom,file=intname,status='unknown',access='append')
3867        open(ipdb,file=pdbname,status='unknown')
3868        open(imol2,file=mol2name,status='unknown')
3869        open(istat,file=statname,status='unknown',access='append')
3870       else
3871 !1out       open(iout,file=outname,status='unknown')
3872       endif
3873 #endif
3874       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3875       csa_seed=prefix(:lenpre)//'.CSA.seed'
3876       csa_history=prefix(:lenpre)//'.CSA.history'
3877       csa_bank=prefix(:lenpre)//'.CSA.bank'
3878       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3879       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3880       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3881 !!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3882       csa_int=prefix(:lenpre)//'.int'
3883       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3884       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3885       csa_in=prefix(:lenpre)//'.CSA.in'
3886 !      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3887 ! Write file names
3888       if (me.eq.king)then
3889       write (iout,'(80(1h-))')
3890       write (iout,'(30x,a)') "FILE ASSIGNMENT"
3891       write (iout,'(80(1h-))')
3892       write (iout,*) "Input file                      : ",&
3893         pref_orig(:ilen(pref_orig))//'.inp'
3894       write (iout,*) "Output file                     : ",&
3895         outname(:ilen(outname))
3896       write (iout,*)
3897       write (iout,*) "Sidechain potential file        : ",&
3898         sidename(:ilen(sidename))
3899 #ifndef OLDSCP
3900       write (iout,*) "SCp potential file              : ",&
3901         scpname(:ilen(scpname))
3902 #endif
3903       write (iout,*) "Electrostatic potential file    : ",&
3904         elename(:ilen(elename))
3905       write (iout,*) "Cumulant coefficient file       : ",&
3906         fouriername(:ilen(fouriername))
3907       write (iout,*) "Torsional parameter file        : ",&
3908         torname(:ilen(torname))
3909       write (iout,*) "Double torsional parameter file : ",&
3910         tordname(:ilen(tordname))
3911       write (iout,*) "SCCOR parameter file : ",&
3912         sccorname(:ilen(sccorname))
3913       write (iout,*) "Bond & inertia constant file    : ",&
3914         bondname(:ilen(bondname))
3915       write (iout,*) "Bending parameter file          : ",&
3916         thetname(:ilen(thetname))
3917       write (iout,*) "Rotamer parameter file          : ",&
3918         rotname(:ilen(rotname))
3919 !el----
3920 #ifndef CRYST_THETA
3921       write (iout,*) "Thetpdb parameter file          : ",&
3922         thetname_pdb(:ilen(thetname_pdb))
3923 #endif
3924 !el
3925       write (iout,*) "Threading database              : ",&
3926         patname(:ilen(patname))
3927       if (lentmp.ne.0) &
3928       write (iout,*)" DIRTMP                          : ",&
3929         tmpdir(:lentmp)
3930       write (iout,'(80(1h-))')
3931       endif
3932       return
3933       end subroutine openunits
3934 !-----------------------------------------------------------------------------
3935       subroutine readrst
3936
3937       use geometry_data, only: nres,dc
3938       use energy_data, only: usampl,iset
3939       use MD_data
3940 !      implicit real*8 (a-h,o-z)
3941 !      include 'DIMENSIONS'
3942 !      include 'COMMON.CHAIN'
3943 !      include 'COMMON.IOUNITS'
3944 !      include 'COMMON.MD'
3945 !el local variables
3946       integer ::i,j
3947 !     real(kind=8) :: var,ene
3948
3949       open(irest2,file=rest2name,status='unknown')
3950       read(irest2,*) totT,EK,potE,totE,t_bath
3951       do i=1,2*nres
3952          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
3953       enddo
3954       do i=1,2*nres
3955          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
3956       enddo
3957       if(usampl) then
3958              read (irest2,*) iset
3959       endif
3960       close(irest2)
3961       return
3962       end subroutine readrst
3963 !-----------------------------------------------------------------------------
3964       subroutine read_fragments
3965
3966       use energy_data
3967 !      use geometry
3968       use control_data, only:out1file
3969       use MD_data
3970       use MPI_data
3971 !      implicit real*8 (a-h,o-z)
3972 !      include 'DIMENSIONS'
3973 #ifdef MPI
3974       include 'mpif.h'
3975 #endif
3976 !      include 'COMMON.SETUP'
3977 !      include 'COMMON.CHAIN'
3978 !      include 'COMMON.IOUNITS'
3979 !      include 'COMMON.MD'
3980 !      include 'COMMON.CONTROL'
3981 !el local variables
3982       integer :: i
3983 !     real(kind=8) :: var,ene
3984
3985       read(inp,*) nset,nfrag,npair,nfrag_back
3986
3987 !el from module energy
3988 !      if(.not.allocated(mset)) allocate(mset(nset))  !(maxprocs/20)
3989       if(.not.allocated(wfrag_back)) then
3990         allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3991         allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
3992
3993         allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
3994         allocate(ifrag(2,nfrag,nset))  !(2,50,maxprocs/20)
3995
3996         allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
3997         allocate(ipair(2,npair,nset))  !(2,100,maxprocs/20)
3998       endif
3999
4000       if(me.eq.king.or..not.out1file) &
4001        write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4002         " nfrag_back",nfrag_back
4003       do iset=1,nset
4004          read(inp,*) mset(iset)
4005        do i=1,nfrag
4006          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4007            qinfrag(i,iset)
4008          if(me.eq.king.or..not.out1file) &
4009           write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4010            ifrag(2,i,iset), qinfrag(i,iset)
4011        enddo
4012        do i=1,npair
4013         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4014           qinpair(i,iset)
4015         if(me.eq.king.or..not.out1file) &
4016          write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4017           ipair(2,i,iset), qinpair(i,iset)
4018        enddo
4019        do i=1,nfrag_back
4020         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4021            wfrag_back(3,i,iset),&
4022            ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4023         if(me.eq.king.or..not.out1file) &
4024          write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4025          wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4026        enddo
4027       enddo
4028       return
4029       end subroutine read_fragments
4030 !-----------------------------------------------------------------------------
4031 ! shift.F       io_csa
4032 !-----------------------------------------------------------------------------
4033       subroutine csa_read
4034
4035       use csa_data
4036 !      implicit real*8 (a-h,o-z)
4037 !      include 'DIMENSIONS'
4038 !      include 'COMMON.CSA'
4039 !      include 'COMMON.BANK'
4040 !      include 'COMMON.IOUNITS'
4041 !el local variables
4042 !     integer :: ntf,ik,iw_pdb
4043 !     real(kind=8) :: var,ene
4044
4045       open(icsa_in,file=csa_in,status="old",err=100)
4046        read(icsa_in,*) nconf
4047        read(icsa_in,*) jstart,jend
4048        read(icsa_in,*) nstmax
4049        read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4050        read(icsa_in,*) nran0,nran1,irr
4051        read(icsa_in,*) nseed
4052        read(icsa_in,*) ntotal,cut1,cut2
4053        read(icsa_in,*) estop
4054        read(icsa_in,*) icmax,irestart
4055        read(icsa_in,*) ntbankm,dele,difcut
4056        read(icsa_in,*) iref,rmscut,pnccut
4057        read(icsa_in,*) ndiff
4058       close(icsa_in)
4059
4060       return
4061
4062  100  continue
4063       return
4064       end subroutine csa_read
4065 !-----------------------------------------------------------------------------
4066       subroutine initial_write
4067
4068       use csa_data
4069 !      implicit real*8 (a-h,o-z)
4070 !      include 'DIMENSIONS'
4071 !      include 'COMMON.CSA'
4072 !      include 'COMMON.BANK'
4073 !      include 'COMMON.IOUNITS'
4074 !el local variables
4075 !     integer :: ntf,ik,iw_pdb
4076 !     real(kind=8) :: var,ene
4077
4078       open(icsa_seed,file=csa_seed,status="unknown")
4079        write(icsa_seed,*) "seed"
4080       close(31)
4081 #if defined(AIX) || defined(PGI)
4082        open(icsa_history,file=csa_history,status="unknown",&
4083         position="append")
4084 #else
4085        open(icsa_history,file=csa_history,status="unknown",&
4086         access="append")
4087 #endif
4088        write(icsa_history,*) nconf
4089        write(icsa_history,*) jstart,jend
4090        write(icsa_history,*) nstmax
4091        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4092        write(icsa_history,*) nran0,nran1,irr
4093        write(icsa_history,*) nseed
4094        write(icsa_history,*) ntotal,cut1,cut2
4095        write(icsa_history,*) estop
4096        write(icsa_history,*) icmax,irestart
4097        write(icsa_history,*) ntbankm,dele,difcut
4098        write(icsa_history,*) iref,rmscut,pnccut
4099        write(icsa_history,*) ndiff
4100
4101        write(icsa_history,*)
4102       close(icsa_history)
4103
4104       open(icsa_bank1,file=csa_bank1,status="unknown")
4105        write(icsa_bank1,*) 0
4106       close(icsa_bank1)
4107
4108       return
4109       end subroutine initial_write
4110 !-----------------------------------------------------------------------------
4111       subroutine restart_write
4112
4113       use csa_data
4114 !      implicit real*8 (a-h,o-z)
4115 !      include 'DIMENSIONS'
4116 !      include 'COMMON.IOUNITS'
4117 !      include 'COMMON.CSA'
4118 !      include 'COMMON.BANK'
4119 !el local variables
4120 !     integer :: ntf,ik,iw_pdb
4121 !     real(kind=8) :: var,ene
4122
4123 #if defined(AIX) || defined(PGI)
4124        open(icsa_history,file=csa_history,position="append")
4125 #else
4126        open(icsa_history,file=csa_history,access="append")
4127 #endif
4128        write(icsa_history,*)
4129        write(icsa_history,*) "This is restart"
4130        write(icsa_history,*)
4131        write(icsa_history,*) nconf
4132        write(icsa_history,*) jstart,jend
4133        write(icsa_history,*) nstmax
4134        write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4135        write(icsa_history,*) nran0,nran1,irr
4136        write(icsa_history,*) nseed
4137        write(icsa_history,*) ntotal,cut1,cut2
4138        write(icsa_history,*) estop
4139        write(icsa_history,*) icmax,irestart
4140        write(icsa_history,*) ntbankm,dele,difcut
4141        write(icsa_history,*) iref,rmscut,pnccut
4142        write(icsa_history,*) ndiff
4143        write(icsa_history,*)
4144        write(icsa_history,*) "irestart is: ", irestart
4145
4146        write(icsa_history,*)
4147       close(icsa_history)
4148
4149       return
4150       end subroutine restart_write
4151 !-----------------------------------------------------------------------------
4152 ! test.F
4153 !-----------------------------------------------------------------------------
4154       subroutine write_pdb(npdb,titelloc,ee)
4155
4156 !      implicit real*8 (a-h,o-z)
4157 !      include 'DIMENSIONS'
4158 !      include 'COMMON.IOUNITS'
4159       character(len=50) :: titelloc1
4160       character*(*) :: titelloc
4161       character(len=3) :: zahl
4162       character(len=5) :: liczba5
4163       real(kind=8) :: ee
4164       integer :: npdb   !,ilen
4165 !el      external ilen
4166 !el local variables
4167       integer :: lenpre
4168 !     real(kind=8) :: var,ene
4169
4170       titelloc1=titelloc
4171       lenpre=ilen(prefix)
4172       if (npdb.lt.1000) then
4173        call numstr(npdb,zahl)
4174        open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4175       else
4176         if (npdb.lt.10000) then
4177          write(liczba5,'(i1,i4)') 0,npdb
4178         else
4179          write(liczba5,'(i5)') npdb
4180         endif
4181         open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4182       endif
4183       call pdbout(ee,titelloc1,ipdb)
4184       close(ipdb)
4185       return
4186       end subroutine write_pdb
4187 !-----------------------------------------------------------------------------
4188 ! thread.F
4189 !-----------------------------------------------------------------------------
4190       subroutine write_thread_summary
4191 ! Thread the sequence through a database of known structures
4192       use control_data, only: refstr
4193 !      use geometry
4194       use energy_data, only: n_ene_comp
4195       use compare_data
4196 !      implicit real*8 (a-h,o-z)
4197 !      include 'DIMENSIONS'
4198 #ifdef MPI
4199       use MPI_data      !include 'COMMON.INFO'
4200       include 'mpif.h'
4201 #endif
4202 !      include 'COMMON.CONTROL'
4203 !      include 'COMMON.CHAIN'
4204 !      include 'COMMON.DBASE'
4205 !      include 'COMMON.INTERACT'
4206 !      include 'COMMON.VAR'
4207 !      include 'COMMON.THREAD'
4208 !      include 'COMMON.FFIELD'
4209 !      include 'COMMON.SBRIDGE'
4210 !      include 'COMMON.HEADER'
4211 !      include 'COMMON.NAMES'
4212 !      include 'COMMON.IOUNITS'
4213 !      include 'COMMON.TIME1'
4214
4215       integer,dimension(maxthread) :: ip
4216       real(kind=8),dimension(0:n_ene) :: energia
4217 !el local variables
4218       integer :: i,j,ii,jj,ipj,ik,kk,ist
4219       real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4220
4221       write (iout,'(30x,a/)') &
4222        '  *********** Summary threading statistics ************'
4223       write (iout,'(a)') 'Initial energies:'
4224       write (iout,'(a4,2x,a12,14a14,3a8)') &
4225        'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4226        'RMSnat','NatCONT','NNCONT','RMS'
4227 ! Energy sort patterns
4228       do i=1,nthread
4229         ip(i)=i
4230       enddo
4231       do i=1,nthread-1
4232         enet=ener(n_ene-1,ip(i))
4233         jj=i
4234         do j=i+1,nthread
4235           if (ener(n_ene-1,ip(j)).lt.enet) then
4236             jj=j
4237             enet=ener(n_ene-1,ip(j))
4238           endif
4239         enddo
4240         if (jj.ne.i) then
4241           ipj=ip(jj)
4242           ip(jj)=ip(i)
4243           ip(i)=ipj
4244         endif
4245       enddo
4246       do ik=1,nthread
4247         i=ip(ik)
4248         ii=ipatt(1,i)
4249         ist=nres_base(2,ii)+ipatt(2,i)
4250         do kk=1,n_ene_comp
4251           energia(i)=ener0(kk,i)
4252         enddo
4253         etot=ener0(n_ene_comp+1,i)
4254         rmsnat=ener0(n_ene_comp+2,i)
4255         rms=ener0(n_ene_comp+3,i)
4256         frac=ener0(n_ene_comp+4,i)
4257         frac_nn=ener0(n_ene_comp+5,i)
4258
4259         if (refstr) then
4260         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4261         i,str_nam(ii),ist+1,&
4262         (energia(print_order(kk)),kk=1,nprint_ene),&
4263         etot,rmsnat,frac,frac_nn,rms
4264         else
4265         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4266         i,str_nam(ii),ist+1,&
4267         (energia(print_order(kk)),kk=1,nprint_ene),etot
4268         endif
4269       enddo
4270       write (iout,'(//a)') 'Final energies:'
4271       write (iout,'(a4,2x,a12,17a14,3a8)') &
4272        'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4273        'RMSnat','NatCONT','NNCONT','RMS'
4274       do ik=1,nthread
4275         i=ip(ik)
4276         ii=ipatt(1,i)
4277         ist=nres_base(2,ii)+ipatt(2,i)
4278         do kk=1,n_ene_comp
4279           energia(kk)=ener(kk,ik)
4280         enddo
4281         etot=ener(n_ene_comp+1,i)
4282         rmsnat=ener(n_ene_comp+2,i)
4283         rms=ener(n_ene_comp+3,i)
4284         frac=ener(n_ene_comp+4,i)
4285         frac_nn=ener(n_ene_comp+5,i)
4286         write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4287         i,str_nam(ii),ist+1,&
4288         (energia(print_order(kk)),kk=1,nprint_ene),&
4289         etot,rmsnat,frac,frac_nn,rms
4290       enddo
4291       write (iout,'(/a/)') 'IEXAM array:'
4292       write (iout,'(i5)') nexcl
4293       do i=1,nexcl
4294         write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4295       enddo
4296       write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4297        'Max. time for threading step ',max_time_for_thread,&
4298        'Average time for threading step: ',ave_time_for_thread
4299       return
4300       end subroutine write_thread_summary
4301 !-----------------------------------------------------------------------------
4302       subroutine write_stat_thread(ithread,ipattern,ist)
4303
4304       use energy_data, only: n_ene_comp
4305       use compare_data
4306 !      implicit real*8 (a-h,o-z)
4307 !      include "DIMENSIONS"
4308 !      include "COMMON.CONTROL"
4309 !      include "COMMON.IOUNITS"
4310 !      include "COMMON.THREAD"
4311 !      include "COMMON.FFIELD"
4312 !      include "COMMON.DBASE"
4313 !      include "COMMON.NAMES"
4314       real(kind=8),dimension(0:n_ene) :: energia
4315 !el local variables
4316       integer :: ithread,ipattern,ist,i
4317       real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4318
4319 #if defined(AIX) || defined(PGI)
4320       open(istat,file=statname,position='append')
4321 #else
4322       open(istat,file=statname,access='append')
4323 #endif
4324       do i=1,n_ene_comp
4325         energia(i)=ener(i,ithread)
4326       enddo
4327       etot=ener(n_ene_comp+1,ithread)
4328       rmsnat=ener(n_ene_comp+2,ithread)
4329       rms=ener(n_ene_comp+3,ithread)
4330       frac=ener(n_ene_comp+4,ithread)
4331       frac_nn=ener(n_ene_comp+5,ithread)
4332       write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4333         ithread,str_nam(ipattern),ist+1,&
4334         (energia(print_order(i)),i=1,nprint_ene),&
4335         etot,rmsnat,frac,frac_nn,rms
4336       close (istat)
4337       return
4338       end subroutine write_stat_thread
4339 !-----------------------------------------------------------------------------
4340 #endif
4341 !-----------------------------------------------------------------------------
4342       end module io_config