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