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