source/unres/io_config.f90

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